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Preface 


This report presents the latest in a number of versions of chemical 
equilibrium and applications programs developed at the NASA Lewis Research 
Center over more than 40 years. These programs have changed over the years to 
include additional features and improved calculation techniques and to take 
advantage of constantly improving computer capabilities. The minimization-of- 
free-energy approach to chemical equilibrium calculations has been used in all 
versions of the program since 1967. 

The two principal purposes of this report are presented in two parts. The 
first purpose, which is accomplished here in part I, is to present in detail a 
number of topics of general interest in complex equilibrium calculations. These 
topics include mathematical analyses and techniques for obtaining chemical 
equilibrium; formulas for obtaining thermodynamic and transport mixture 
properties and thermodynamic derivatives; criteria for inclusion of condensed 
phases; calculations at a triple point; inclusion of ionized species; and various 
applications, such as constant-pressure or constant-volume combustion, rocket 
performance based on either a finite- or infinite-chamber-area model, shock wave 
calculations, and Chapman-Jouguet detonations. 

The second purpose of this report, to facilitate the use of the computer 
code, is accomplished in part II, entitled “Users Manual and Program 
Description.” Various aspects of the computer code are discussed, and a number 
of examples are given to illustrate its versatility. 


p«ao**N* page plank not filmed 

/ 

iii 


PAGE— Ll— INTENTIONALLY 8Ur< 





Contents 

Page 

1. Introduction 1 

2. Equations Describing Chemical Equilibrium 3 

2.1 Units 3 

2.2 Equation of State 3 

2.3 Minimization of Gibbs Energy 4 

2.3.1 Gibbs Iteration Equations 6 

2.3.2 Reduced Gibbs Iteration Equations 7 

2.4 Minimization of Helmholtz Energy 8 

2.4.1 Helmholtz Iteration Equations 9 

2.4.2 Reduced Helmholtz Iteration Equations 9 

2.5 Thermodynamic Derivatives From Matrix Solutions 10 

2.5.1 Derivatives With Respect to Temperature 10 

2.5.2 Derivatives With Respect to Pressure 11 

2.6 Other Thermodynamic Derivatives 12 

3. Procedure for Obtaining Equilibrium Compositions 13 

3.1 Initial Estimates 13 

3.2 Magnitude of Species Used During Iteration 13 

3.3 Convergence 14 

3.4 Tests for Condensed Phases 15 

3.5 Phase Transitions and Special Derivatives 16 

3.6 Singularities 17 

3.7 Iteration Procedure and Tests for Ions 17 

4. Thermodynamic Data 19 

4.1 Assigned Enthalpies 19 

4.2 Least-Squares Coefficients 19 

5. Thermal Transport Property Data 21 

5. 1 Data for Individual Species 21 

5.2 Mixture Property Data 21 

5.2.1 Viscosity and Frozen Thermal Conductivity 21 

5.2.2 Reaction Thermal Conductivity 22 

5.2.3 Specific Heat for Gases Only 23 

5.2.4 Prandtl Number 23 

6. Theoretical Rocket Performance 25 

6.1 Assumptions 25 

pnftOi&ttMl PAGE tLANfc NOT FILMED v p ^\1 


INTENTIONALLY BLANK 


6.2 Parameters 26 

6.2.1 Conservation Equations 26 

6.2.2 Velocity of Flow 26 

6.2.3 Force 26 

6.2.4 Specific Impulse 26 

6.2.5 Mach Number 27 

6.2.6 Characteristic Velocity 27 

6.2.7 Area per Unit Mass Flow Rate 27 

6.2.8 Coefficient of Thrust 27 

6.2.9 Area Ratio 27 

6.3 Procedure for Obtaining Equilibrium Rocket Performance 

for IAC Model 27 

6.3.1 Combustion Conditions 28 

6.3.2 Exit Conditions 28 

6.3.3 Throat Conditions 28 

6.3.4 Discontinuities at Throat 28 

6.3.5 Assigned Subsonic or Supersonic Area Ratios 29 

6.3.6 Empirical Formulas for Initial Estimates of P m JP e 29 

6.3.7 Analytic Expression for Improved Estimates 

ofPJP, 29 

6.4 Procedure for Obtaining Equilibrium Rocket Performance 

for FAC Model 30 

6.4. 1 Initial Estimates for P mi and AJA t 30 

6.4.2 Improved Estimates for ^inf and AJA t 31 

6.5 Procedure for Obtaining Frozen Rocket Performance 31 

6.5.1 Exit Conditions 31 

6.5.2 Throat Conditions 32 

6.5.3 Thermodynamic Derivatives for Frozen Composition ... 32 

7. Incident and Reflected Shocks 33 

7.1 Incident Shocks 33 

7.1.1 Iteration Equations 34 

7.1.2 Corrections and Convergence 35 

7.1.3 Initial Estimates of T^T X and P 2 /P\ 35 

7.2 Reflected Shocks 36 

7.2.1 Iteration Equations 36 

7.2.2 Corrections and Convergence 37 

7.2.3 Initial Estimates of T 5 IT 2 and P 5 /P 2 37 

8. Chapman-Jouguet Detonations 39 

8.1 Iteration Equations 39 

8.2 Corrections and Convergence 40 

8.3 Initial Estimates of T 2 IT X and P 2 /P\ 40 

9. Input Calculations 41 

Appendix — Symbols 45 

References 49 


vi 



Chapter 1 

Introduction 


Knowing the chemical equilibrium compositions of 
a chemical system permits one to calculate theoretical 
thermodynamic properties for the system. These proper- 
ties can be applied to a wide variety of problems in 
chemistry and chemical engineering. Some applications 
are the design and analysis of equipment such as 
compressors, turbines, nozzles, engines, shock tubes, heat 
exchangers, and chemical processing equipment. 

For more than 40 years the NASA Lewis Research 
Center has been involved in developing methods and 
computer programs for calculating complex chemical 
equilibrium compositions and thermodynamic properties 
of the equilibrium mixtures and in applying these 
properties to a number of problems (Gordon et al., 1959, 
1962, 1963, 1970, 1976, 1984, 1988; Huff et al., 1951; 
Svehla and McBride, 1973; and Zeleznik and Gordon, 
1960, 1962a, b, 1968). Earlier versions of the chemical 
equilibrium computer program (Zeleznik and Gordon, 
1962a; and Gordon and McBride, 1976) have had wide 
acceptance. Since the last publication this program has 
been under continuous revision and updating, including 
several substantial additions. One addition is an option 
for obtaining the transport properties of complex 
mixtures (Gordon et al., 1984) by methods simpler than 
those of Svehla and McBride (1973). A second addition 
is an option to calculate rocket performance for a rocket 
motor with a finite-area combustor (Gordon and 
McBride, 1988). The present report documents these and 
other additions and revisions to the program since 1976. 
The revised program is called CEA (Chemical 
Equilibrium and Applications). 

The program can now do the following kinds of 
problems: 

(1) Obtaining chemical equilibrium compositions 
for assigned thermodynamic states. These states may be 
specified by assigning two thermodynamic state func- 
tions as follows: 


(a) Temperature and pressure, tp 

(b) Enthalpy and pressure, hp 

(c) Entropy and pressure, sp 

(d) Temperature and volume (or density), tv 

(e) Internal energy and volume (or density), uv 

(f) Entropy and volume (or density), sv 

(2) Calculating theoretical rocket performance for 
a finite- or infinite-area combustion chamber 

(3) Calculating Chapman-Jouguet detonations 

(4) Calculating shock tube parameters for both inci- 
dent and reflected shocks 

Some problems handled by the program use just one 
combination of assigned states — namely, the tp, hp, sp, 
tv, uv, and sv problems. For example, the tp problem, 
which consists of a schedule of one or more assigned 
temperatures and pressures, might be used to construct 
Mollier diagrams. The hp problem gives constant- 
pressure combustion properties and the uv problem gives 
constant-volume combustion properties. Other problems 
make use of more than one combination of assigned 
thermodynamic states. For example, the shock and 
detonation problems use hp and tp; the rocket problem 
uses hp or tp and also sp. 

This report consists of two parts. Part I, containing 
the analysis, includes 

(1) The equations describing chemical equilibrium 
and the applications previously mentioned (i.e., rocket 
performance, shocks, and Chapman-Jouguet detonations) 

(2) The reduction of these equations to forms 
suitable for mathematical solution by means of iterative 
procedures 

(3) Equations for obtaining thermodynamic and 
transport properties of mixtures 

Part II, a program description and users manual, 
discusses the modular form of the program and briefly 
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describes each subroutine. It also discusses the prepara- 
tion of input, various permitted options, output tables, 
and error messages. A number of examples are also 
given to illustrate the versatility of the program. 

In addition to the work in chemical equilibrium 
calculations and applications over the past 40 years, 
progress in computer programs, data generation, and data 
fitting has also been made at NASA Lewis for the ther- 
modynamic and thermal transport properties of indi- 
vidual species required for the equilibrium calculations. 
Some examples of this effort are Burcat et al. (1985), 
McBride et al. (1961, 1963, 1967, 1992, 1993a, b), 


Svehla (1962), and Zeleznik and Gordon (1961). In 
addition to data calculated by us, other thermodynamic 
data included in our files are taken from sources such as 
Chase et al. (1985), Garvin et al. (1987), Gurvich et al. 
(1989), and Marsh et al. (1988). Files distributed with 
the computer program are described in part II. 

Various versions of the equilibrium program or 
modifications of the program have been incorporated 
into a number of other computer codes. An example 
is Radhakrishnan and Bittker (1994) for kinetics cal- 
culations. 
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Chapter 2 

Equations Describing Chemical Equilibrium 


Chemical equilibrium is usually described by either 
of two equivalent formulations — equilibrium constants or 
minimization of free energy. Reports by Zeleznik and 
Gordon (1960, 1968) compare the two formulations. 
Zeleznik and Gordon (1960) shows that, if a generalized 
method of solution is used, the two formulations reduce 
to the same number of iteration equations. However, 
with the minimization-of-free-energy method each 
species can be treated independently without specifying 
a set of reactions a priori, as is required with equilibrium 
constants. Therefore, the minimization-of-free-energy 
formulation is used in the CEA program. 

The condition for equilibrium can be stated in terms 
of any of several thermodynamic functions, such as the 
minimization of Gibbs or Helmholtz energy or the 
maximization of entropy. If one wishes to use tempera- 
ture and pressure to characterize a thermodynamic state, 
Gibbs energy is most easily minimized inasmuch as 
temperature and pressure are its natural variables. 
Similarly, Helmholtz energy is most easily minimized if 
the thermodynamic state is characterized by temperature 
and volume (or density). 

Zeleznik and Gordon (1960) presents equations 
based on minimization of Gibbs energy. Some of these 
equations are repeated and expanded herein for 
convenience. In addition, a set of equations based on 
minimization of Helmholtz energy is also presented. 
However, because only ideal gases and pure condensed 
phases are being considered, the general notation of 
Zeleznik and Gordon (1960) is not used. 

2.1 Units 

The International System of Units (SI) used in this 
report is 


Physical 

quantity 

Unit 

Symbol 

Length 

meter 

m 

Mass 

kilogram 

kg 

Time 

second 

s 

Temperature 

kelvin 

K 

Force 

newton 

N 

Pressure 

newton per 
square meter 

N/m 2 

Work, energy 

joule 

J 


The numerical values of a number of fundamental 
contants are taken from Cohen and Taylor (1987). For 
example, the value of the gas constant R taken from this 
reference is 8314.51 J/(kg-mole)(K). In those sections 
dealing with the computer program, other units are used 
in addition to or instead of SI units. 

2.2 Equation of State 

In this report we assume that all gases are ideal and 
that interactions among phases may be neglected. The 
equation of state for the mixture is 

PV = nRT ( 2 - la ) 

or 

l = nRT (2-lb) 

P 
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where P is pressure (in newtons per square meter), V 
specific volume (in cubic meters pier kilogram), n moles 
per unit mass of mixture (in kilograms-mole per 
kilogram), T temperature (in kelvin), and p density (in 
kilograms per cubic meter). Symbols used in this report 
are defined in the appendix. For a reacting chemical 
system the number of moles n is generally not constant. 

Equation (2.1) is assumed to be correct even when 
small amounts of condensed species (up to several 
percent by weight) are present. In this event the 
condensed species are assumed to occupy a negligible 
volume relative to the gaseous species. An example 
given in part II of this report illustrates the validity of 
this assumption. In the variables V, n, and p the volume 
and mole number refer to gases only, but the mass is for 
the entire mixture including condensed species. The 
word “mixture” is used in this report to refer to mixtures 
of reaction products as distinguished from mixtures of 
reactants, which are referred to as “total reactants.” 

On the basis of this definition, n can be written as 

NG 

n = ][>, (2-2) 

/-i 


MW = 


M 



(2.4a) 


Molecular weight is given the symbol MW in equa- 
tion (2.4a) to differentiate it from M, The two different 
definitions of molecular weight, M and MW, give differ- 
ent results only in mixtures of products containing con- 
densed as well as gaseous species. Only M is given in 
the output, but MW may be obtained from M by means 
of 


MW =M 


NS 


i- E 


(2.4b) 


/=NG+1 J 


where x- is the mole fraction of species j relative to all 
species in the mixture. Some additional discussion of the 
differences in these molecular weights is given in part II 
of this report. 


where rij is the number of kilogram-moles of species j 
per kilogram of mixture and the index NG refers to the 
number of gases in the mixture. The molecular weight of 
the mixture M is defined as 


2.3 Minimization of Gibbs Energy 

For a mixture of NS species the Gibbs energy per 
kilogram of mixture g is given by 


M = — (2.3a) 

n 


NS 

g = E m 


(2.5) 


or equivalently as 


M = 


E"yW, 

tl 

NG 

/=! 


(2.3b) 


where M is the molecular weight of species j and the 
index NS refers to the number of species in the mixture. 
In the CEA computer program, among the NS species, 
gases are indexed from 1 to NG and condensed species 
from NG + 1 to NS. 

More conventionally, molecular weight is defined 
as 


where the chemical potential per kilogram-mole of 
species j is defined to be 


= 


dg_ 

3/i. j 


( 2 . 6 ) 


The condition for chemical equilibrium is the minimiza- 
tion of free energy. This minimization is usually subject 
to certain constraints, such as the following mass-balance 
constraints: 


NS 

Ev, 




- b° 


= 0 


(i = !,...,() (2.7a) 
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or 


b , - b° =0 (i = (2.7b) 

where the stoichiometric coefficients a tJ are the number 
of kilogram-atoms of element i per kilogram-mole of 
species j, the index f is the number of chemical ele- 
ments (if ions are considered, the number of chemical 
elements plus one), b° is the assigned number of 
kilogram-atoms of element i per kilogram of total 
reactants (see eq. (9.5)), and 

NS 

b i = Y, a n n j (‘ = 1 ’ -*0 (2.7c) 

P i 

is the number of kilogram-atoms of element i per kilo- 
gram of mixture. 

Defining a term G to be 

G m 8 + E H b i ~ *?) ( 2 . 8 ) 

i-1 


where X- are Lagrangian multipliers, the condition for 
equilibrium becomes 


NS 


« G = E| n, + E Vvl *"/ + E(*« - - 0 

(2.9) 


J = 1 


i= 1 




Treating the variations 5 rij and 5A, f as independent gives 


»j + E K a ij = 0 


(j * 1.-J9S) (2.10) 


and also gives the mass-balance equation (2.7b). 

From the assumptions in section 2.2 the chemical 
potential can be written as 


state. For a gas the standard state is the hypothetical 
ideal gas at the standard-state pressure. For a pure solid 
or liquid the standard state is the substance in the 
condensed phase at the standard-state pressure. 
Historically, the defined standard-state pressure has been 
1 atmosphere (101 325 Pa). Most early tabulations of 
thermodynamic data were based on this pressure. 
However, in 1982 the International Union of Pure and 
and Applied Chemistry (Cox, 1982) recommended that 
the standard-state pressure should be defined as 1 bar 
(10 5 Pa). Most recent compilations have used 1 bar as 
the standard pressure (e.g., Chase et al. (1985), Garvin 
et al. (1987), Gurvich et al. (1989), Marsh et al. (1988), 
and McBride et al. (1993a, b)). The unit of pressure in 
equation (2.11) should be consistent with the unit of 
pressure in the thermodynamic data being used. 

The term p° and other thermodynamic terms that 
appear later in the text, such as C°j, HJ, Sj, C°j , and 
UJ, are all functions of temperature. However, including 
T as part of the symbol notation, such as H°AT) or H°(T), 
is done only when needed for clarity. For example, 
sensible heat for species j between a temperature T and 
a temperature of 298.15 K can be written as H°(T) - 
tfj(298.l5). 

Equations (2.7a) and (2.10) permit the determin- 
ation of equilibrium compositions for thermodynamic 
states specified by an assigned temperature T Q and 
pressure P Q . That is, in addition to equations (2.7a) and 
(2.10), we have the pair of trivial equations 


II 

(2.12a) 

p = p 
r r o 

(2.12b) 


However, the thermodynamic state can be specified by 
assigning any two state functions. For example, the 
thermodynamic state corresponding to a constant- 
pressure combustion is specified, instead of by 
equations (2.12), by 


= 


\ 1 ° + RT ]n^ + RT In P 
J n 

(j = NG + 1,...,NS) 


,NG) 

II 

*45 

(2.13a) 

(2.11) 

“TJ 

ll 

(2.13b) 


where h is the specific enthalpy of the mixture and h Q a 
where \x° for gases (j = 1 to NG) and for condensed constant equal to the specific enthalpy of the reactants 

phases (j > NG) is the chemical potential in the standard (see eq. (9.7)). The expression for h is 
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NS 

h = E V*/ (214) 

h i 

where /f° is the standard-state molar enthalpy for 
species j at temperature T. 

For assigned entropy and pressure (such as for an 
isentropic compression or expansion to a specified 
pressure), the thermodynamic state is specified by 

s = (2.15a) 


pansion of the appropriate equations with all terms 
truncated that contain derivatives higher than the first. 
The correction variables used are Ain n. (j = 1,..,,NG), 
A rij (j = NG + 1,...,NS), Aln/i, n f = -XJRT, and 

AlnT. As Zeleznik and Gordon (1968) points out, it is 
no restriction to start each iteration with the estimate for 
the Lagrangian multipliers equal to zero inasmuch as 
they appear linearly in equation (2.10). After making 
dimensionless those equations containing thermodynamic 
functions, the Newton-Raphson equations obtained from 
equations (2.10), (2.7), (2.13a), and (2.15a) are 


P = P 0 (2.15b) 


where s is the specific entropy of the mixture and s 0 the 
assigned specific entropy, or the specific entropy of the 
total reactant (see eq. (9.22)). The expression for s is 


s 


NS 


= E"/ 5 , 


(2.16) 


where 


H? 


Ain n, - Y - Aln/i — Ain 7 

J £—1 V ' DT 


i=l 


RT 




= U = 1,. .^G) 


(2.18) 


-Ew 


i=l 


H° 

— Ain 7 = 
RT 


RT 


O’ = NG + 1,...,NS) (2.19) 



- R In ^ - R]n P 

n 

(j = NG + 1.....NS) 


U 


1,..,NG) 

(2.17) 


and S° is the standard- state molar entropy for species j. 
Equation (2.17) is similar to equation (2.11), and the 
same discussion concerning standard-state pressure that 
applied to equation (2. 1 1) also applies to equation (2. 17). 

2.3.1 Gibbs Iteration Equations 

The equations required to obtain composition are 
not all linear in the composition variables and therefore 
an iteration procedure is generally required. In the itera- 
tion procedure to be described it will be convenient to 
treat n as an independent variable. A descent Newton- 
Raphson method is used to solve for corrections to initial 
estimates of compositions n , Lagrangian multipliers A. , 
moles of gaseous species n, and (when required) 
temperature 7. This method involves a Taylor series ex- 


NG NS 

EVVAlnV E % An j = b °k ~ b k 

j-\ /= NG + 1 


(k = U.,9 


(2.20) 


NG 


NG 


y ftj Ain rij - n Alnn = n - ^ n . 

>■ l M 


( 2 . 21 ) 



NS nC° 
\j - 1 « 


Ain 7> 


K - b 

RT 


( 2 . 22 ) 
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NG 

E 


/? 


Ain n. 


NS 

► E 

/-NG • 1 


S.° 


' NS 


A n, 


n.C, 


E J 


P.i 


Ain 7 


V* 1 


i-1 


g 

V>*1 


K7 


NS 


*i + 


_ 

E —An. 

^ K7 1 


y-NG*l 


NO 

E 

V«l 


*7 j 


Alnn 


s - s NG 

--V- + *-E», < 2 - 23 ) 

A j- 1 

where C°j is the standard-state specific heat at constant 
pressure for species j at temperature T. 

2.3.2 Reduced Gibbs Iteration Equations 




P R U R 2 T 2 

RT P r 2 t 2 


(2.27) 


For problems with assigned thermodynamic states 
tp, hp , or sp, various combinations of equations (2.18) to 
(2.23) could be used to obtain corrections to estimates. 
However, for chemical systems containing many species, 
it would be necessary to solve a large number of 
simultaneous equations. This large number of equations 
can be reduced quite simply to a much smaller number 
by algebraic substitution. The expression for Ain n. 
obtained from equation (2.18) is substituted into equa- 
tions (2.20) to (2.23). When equation (2.19) written with 
signs reversed is included, the resulting reduced 
equations are 


I NG 

EEw/i 


i=i j-i 


NS /NG 

E a n An j + E<v, 

/=NG+1 V= 1 


Ain n 


f NG n w 

£ V/v 


NG 


V-i 


RT 


Ain r = + 


(Jt - 


(2.24) 


H° 


i= 1 


Foi + — Ain 7 = — (j - NG + 1,...,NS) 
^ y ' RT RT 


(2.25) 


I NG 


/NG 


EE w. + E w /-« 

i=l /'=! (/=1 ) 


Ain n 


NG 


V;=i 


/?7 


Ain T 


NG NG 

= »-E«, * E 


*7 


(2.26) 


I NG 

E E 

i-i \/-i 



NS 


♦ E 

>=NG*1 



NS 


0=1 


'7 


n,c:, ™ HjHJSj 

h R 2 T ) 

NG 


J n -5 

Ain 7 = -2 — + n 


5° n.S, p y 

'E n / + E4^ (2-28) 




T - T * z 7 


Equations (2.24) to (2.28) are given in table 2.1 in 
a form that permits a direct comparison with other sets 
of simultaneous equations presented in later sections 
(tables 2.2 to 2.4). (Note that tables 2.1 to 2.4 appear at 
the end of this document.) In a previous report (Gordon 
and McBride, 1976), some special script symbols were 
used for the sake of compactness in preparing tables 2. 1 
to 2.4, as for example, 


Of = — (2-29) 

RT 

These script symbols are no longer used in this report. 

The correction equations required for several types 
of constant-pressure problems are summarized as 
follows, where i = and j = NG + 1,...,NS: 


Type of 
problem 

Equations 

required 

Correction 

variables 

Assigned temperature 
and pressure, tp 

(2.24), (2.25), 
(2.26) 

7l ( , A rij, Ain n 

Assigned enthalpy 

(2.24), (2.25), 

71-, A rij, Ain n , 

and pressure, hp 

(2.26), (2.27) 

Ain T 

Assigned entropy and 

(2.24), (2.25), 

7l ( , A rij, Ain n. 

pressure, sp 

(2.26), (2.28) 

Ain T 
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After obtaining the correction variables shown 
above, the corrections for gaseous species Ain rij 

(j = 1.....NG) are then obtained from equation (2.18). 
Section 3.3 discusses controlling the size of corrections 
before they are applied to obtain improved estimates. 


2.4 Minimization of Helmholtz Energy 


NS 

y-i 


+ Y,{ b i ~ b i°) iX i = 0 

i=l ) <*1 

(2.34) 


i 

E 

1*1 


Treating bn- and as independent again gives, as in 
section 2.3, equations (2.7) and (2.10). Now, however, 
instead of equation (2.11), 


The equations presented in this section are similar 
to those in section 2.3. Whatever differences appear are 
due to the different forms of the chemical potential 
\ij (j = 1,...,NG). In section 2.3 pressure was one of the 
assigned thermodynamic states, and consequently Gibbs 
energy was minimized. In this section volume (or 
density) is one of the assigned thermodynamic states, 
and consequently Helmholtz energy is minimized. 

The two energies (Gibbs and Helmholtz) have the 
following thermodynamic relationship: 

/ = g ~ PV (2.30) 

where / is the Helmholtz energy per kilogram of mix- 
ture. After substituting Gibbs energy g as given by 
equation (2.5), equation (2.30) becomes 

NS 

/ = E W ~ PV (2.31) 

The chemical potential \Xj can be expressed as a thermo- 
dynamic derivative in several ways (Kirkwood and 
Oppenheim, 1961). One way is given by equation (2.6). 
Another expression is 


If 


= 


df) 


. dn,\ 

V J >T, v. 




(2.32) 


|ij +/J7'ln 


rtjR'j 


< j = 1.....NG) 


H, = < 


[Vj 


(j = NG + 1,...,NS) 


(2.35) 


where R' = Rx 10 -5 . 

Equations (2.7) rmd (2.10), with p ; given by 
equation (2.35), permit the determination of equilibrium 
compositions for thermodynamic states specified by an 
assigned temperature 7 0 and volume V 0 ; that is, in 
addition to equations (2.7) and (2.10), we have the pair 
of trivial equations 


T=T 0 (2.36a) 

V = V 0 (2.36b) 

Analogous to equation (2.13) for a constant- 
pressure combustion process, we can set down the 
following conditions for constant-volume combustion: 

u' = u £ (2.37a) 

V = V 0 (2.37b) 

where u' is the specific internal energy of the mixture 
and u'q a constant equal to the specific internal energy of 
the reactants. The expression for u' is 


F =/ + EMv fe .°) (2 - 33) 

M 


NS 

*' - E"X 


(2.38) 


the condition for equilibrium based on the minimization 
of Helmholtz energy subject to mass-balance constraints 
is 


where U° is the standard-state molar internal energy for 
species j. 
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Analogous to equation (2.15), for assigned entropy 
and volume (such as for an isentropic compression or 
expansion to a specified volume), the thermodynamic 
state is specified by 


s = s 0 (2.39a) 

V = V 0 (2.39b) 

Iteration equations are derived in section 2.4.1 that 
permit solution of the composition variables for constant- 
volume problems. 

2.4.1 Helmholtz Iteration Equations 

Correction equations are obtained in a manner 
similar to that described in section 2.3.1. In this case, 
however, the expression for \ij is equation (2.35) rather 
than equation (2.11). Because n does not appear 
explicitly as a variable in equation (2.35), Ain n no 
longer appears as a correction variable. The Newton- 
Raphson equations obtained from equations (2.10), (2.7), 
(2.37a), and (2.39a) are 
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where C°j is the standard-state specific heat at constant 
volume for species j at temperature T. 


2.4.2 Reduced Helmholtz Iteration Equations 

Equations (2.40) to (2.44) may be reduced to a 
much smaller set of working correction equations by 
eliminating Ain rij , obtained from equation (2.40), from 
equations (2.42) to (2.44). When equation (2.41) written 
with the signs reversed is included, the resulting reduced 
set of equations is 
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composition is assumed to reach its equilibrium value 
instantaneously, the derivatives are referred to as 
“equilibrium” derivatives. If, on the other hand, reaction 
times are assumed to be infinitely slow, composition 
remains fixed (frozen) and the derivatives are referred to 
as “frozen.” Special subscripts are used to differentiate 
these different conditions only for c , thermal 
conductivity, and Prandtl number. The equilibrium value 
of c p may be expressed as the the sum of a “frozen” 
contribution and a “reaction” contribution as follows: 


Equations (2.45) to (2.48) are given in table 2.2 in a 
form that permits direct comparison with the iteration 
equations in table 2.1 and the derivative equations in 
tables 2.3 and 2.4. 

The correction equations required for several types 
of constant-volume problems are summarized as follows, 
where i = !,...,{ and j = NG + 1,...,NS: 


Type of 
problem 

Equations 

required 

Correction 

variables 

Assigned temperature 

(2.45), (2.46) 

*1 • &*j 

and volume, tv 



Assigned internal 
energy and volume, 

(2.45), (2.46), 
(2.47) 

Kj, A jtj, Ain T 

uv 



Assigned entropy and 
volume, sv 

(2.45), (2.46), 
(2.48) 

7t ( , A /ij, Ain T 


After obtaining the correction variables the 
corrections for gaseous species Ain n } (/=1,...,NG) are 
then obtained from equation (2.40). (See section 3.3 for 
discussion on controlling size of corrections before 
applying them to obtain improved estimates.) 
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The expressions in equations (2.49b) and (2.49c) were 
obtained by differentiating equation (2.14). (Subscripts 
different from those used in equation (2.49a) are used in 
transport property calculations. See section 5.2.3). 

From equation (2.1) 



(2.50) 

(2.51) 


2.5.1 Derivatives With Respect to Temperature 


2.5 Thermodynamic Derivatives From 
Matrix Solutions 

All thermodynamic first derivatives can be ex- 
pressed in terms of any three independent first deriva- 
tives. The Bridgman tables, as tabulated, for example, in 
Lewis and Randall (1961), express first derivatives in 
terms of (dV/dT)p, (dVjdP) T , and (dh/d T) p = c p . We use 
the logarithmic form of the volume derivatives because 
it gives an indication of the extent of chemical reaction 
occurring among the reaction species. These derivatives 
may have more than one value depending on what is 
assumed occurs to composition in a thermodynamic 
process from one condition to another. If, for example, 


The derivatives of rij and n with respect to temper- 
ature are needed to evaluate equations (2.49c) and (2.50). 
These may be obtained by differentiating equations 
(2.10), (2.7), and (2.3), which gives the following: 
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As in the case of the iteration correction equations 
previously discussed, the derivative equations can be 
reduced to a much smaller number of simultaneous 
equations by eliminating (3 In n ; /3 In T) p , obtained from 
equation (2.52), from equations (2.54) and (2.55). When 
equation (2.53) written with the sign reversed is 
included, the resulting reduced number of temperature 
derivative equations is 
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RT 
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Equations (2.56) to (2.58) are given in table 2.3 in 
a form that simplifies the comparison of derivative equa- 
tions with the iteration equations in tables 2.1 and 2.2. 
The values of the derivatives obtained from equa- 
tions (2.56) to (2.58) could be used in equation (2.52) to 
obtain derivatives for the gaseous species (3 In rij/d In T) p 

(j = 1.....NG) and all of the temperature derivatives 
could then be used to evaluate c from equations (2.49). 
However, there is an alternative and much simpler 
procedure for obtaining c p . Substituting 3 In njd In T) p 
obtained from equation (2.52) into equations (2.49) and 
dividing by R yield 



(2.59) 

In equation (2.59) only the temperature derivatives 
obtained directly from solution of equations (2.56) to 
(2.58) are required. Furthermore, all the coefficients in 
equation (2.59) are exactly the coefficients appearing in 
the reduced-enthalpy equation (2.27). The second-last 
term in equation (2.59) is the frozen contribution to 
specific heat (eq. (2.49b)); the remainder of the terms are 
the reaction contribution (eq. (2.49c)). 

2.5.2 Derivatives With Respect to Pressure 

The derivative (31n n/dln P) T can be obtained in a 
manner similar to that described for obtaining derivatives 
with respect to the temperature. Differentiating equa- 
tions (2.10), (2.7), and (2.3) gives 

'^1 .-i 

31nPj T ft VlnPj T UlnPj r 
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3lnP 


T /’NG * 1 


31nP 


(* = 


r,te| -„[««] =o 

U Vtal>j r lainPj, 


Equations (2.60) to (2.63) can be reduced to a 
smaller set by eliminating (31n njd In F) T , obtained from 
equation (2.60), from equations (2.62) and (2.63). When 
equation (2.61) written with the sign reversed is 
included, the results are 
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Using the symbols 


and 



Equations (2.64) to (2.66) are given in table 2.4 for 
comparison with tables 2.1, 2.2, and 2.3. The only 
derivative obtained from solution of equations (2.64) 
to (2.66) that is used is (d]nnJd]nP) T (see eq. (2.51). 


equation (2.69) may be written as 


y s - - 



(2.70) 


(2.71) 


(2.72) 


(2.73) 


2.6 Other Thermodynamic Derivatives 

As stated previously, all thermodynamic first 
derivatives can be expressed in terms of the three 
thermodynamic first derivatives discussed in the previous 
sections — namely, c py 31n Vjd\nT) p , and (din Vjd\nP) T 
(see Bridgman tables in Lewis and Randall (1961)). 
Velocity of sound a , a frequently used parameter, is 
defined by 


«’ . W <2.67) 

UpJs pUtapJ, p(.3to*7j 

From Bridgman tables 



. / 3 In V \ + PV( d In V ? 

H3inpJ T + riainrj. 
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This may be written as 


Using the equation of state given in equation (2.1), 
we obtain from equation (2.67) the familiar expression 
for velocity of sound 


a = s /nRTy s (2.74) 

Note that the y s defined by equation (2.71) is required 
in equation (2.74) and not the specific heat ratio y 
defined in equation (2.72). 

In section 3.5 an alternative expression is derived 
for y s for the special situation of triple phases, where the 
expressions in equations (2.68), (2.69), and (2.73) are no 
longer valid. 

Gordon and Zeleznik (1962) gives numerous first 
derivative relations that are of interest in rocket 
performance calculations. One of these derivatives, 
which is used in section 6.3.4, is 


r d]nP \ 

K d\nT) 


Jain V 
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(2.69) 


A numerical example involving this derivative is given 
in Zeleznik and Gordon (1968). 
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Chapter 3 

Procedure for Obtaining Equilibrium 
Compositions 


In principle, obtaining equilibrium compositions by 
means of the Newton-Raphson iteration procedure dis- 
cussed in sections 2.3.1, 2.3.2, 2.4.1, and 2.4.2 should 
offer no difficulties. However, a number of practical 
items require detailed attention in order to avoid numeri- 
cal difficulties: initial estimates, tests for condensed 
phases, phase transitions and triple points, convergence, 
accidental singularities, special handling of ions, and 
consideration of trace species. 

3.1 Initial Estimates 

An extremely simple procedure is used in this 
report to assign estimates for composition. For the first 
iteration of the first point in a schedule of points, we 
assign n - 0.1, which is equivalent to an estimate of 10 
for molecular weight. Then the number of kilogram- 
moles of each gaseous species per kilogram of mixture 
is set equal to 0.1/NG, where NG is the number of gase- 
ous species being considered. The number of moles of 
each condensed species is set equal to zero. For hp, sp, 
uv, and sv problems an arbitrary initial estimate of 
T = 3800 K is used by the program unless a different 
estimate is included in the input. 

Admittedly, this simple procedure will often give 
poor initial estimates. However, for a general chemistry 
program, we find this technique preferable to the alterna- 
tive of devising numerous special routines for obtaining 
good estimates for numerous possible chemical systems. 
Furthermore, the estimating technique is used only for 
the first point in any schedule of points. For all points 
after the first the results of a preceding point serve as 
initial estimates. 

Because no attempt is made to obtain good initial 
estimates, the question arises whether convergence can 


be “guaranteed.” This question is discussed in 
section 3.3. 

3.2 Magnitude of Species Used During 
Iteration 

Both the linear and logarithmic composition varia- 
bles are used for gaseous species during the composition 
iteration process. Only the linear variable is used for 
condensed species. Corrections to compositions for 
gases are in the form of logarithmic variables Ain /i ■, and 
therefore the logarithmic values of gaseous compositions 
In rij are continuously updated from iteration to iteration. 
The linear values of the compositions rij are obtained by 
taking the antilogarithm of In rij. However, to save com- 
puter time during iteration, n- are calculated only for 
those species whose mole fractions are greater than a 
certain specified size. In a previous version of the 
program (Gordon and McBride, 1976) this specified size 
had only one value, namely n-!n - 10~ 8 (or In nJn 
- -18.420681). A program variable SIZE was defined as 
SIZE = 18.420681. Thus, antilogarithms of In rij were 
obtained only for gases meeting the following condition: 
In nj/n > -SIZE ( rij/n > 1CT 8 .) For gaseous species not 
passing this test rij were set equal to zero. In addition, 
the maximum number of iterations permitted was 35. 

In the present CEA program several additional 
variables relating to SIZE have been added to handle 
more demanding situations. Two of the present variables 
relating to the mole fraction size for which antilogari- 
thms are obtained are TSIZE for non-ionized species and 
ESIZE for ionized species. TSIZE may be modified for 
any of the following reasons: inclusion of species in the 
calculations with mole fractions smaller than 10 -8 (by 
means of an input parameter TRACE); a change of 
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components; a singular matrix; or the chemical system 
under consideration containing a chemical element that 
differs in magnitude from the largest of the other 
elements by more than 10 5 . The purpose of changing 
TSIZE in the last case is to ensure that not all species 
containing the trace element will be eliminated during 
iteration. To aid in testing for trace elements, a param- 
eter BRATIO is defined to be the ratio of the elements 
with the lowest to highest kilogram-atoms per kilogram 
of mixture. The following, then, are the conditions under 
which several parameters relating to species size are set 
to various values. These conditions are based on exten- 
sive practical experience with obtaining solutions for 
many types of problems: 

(1) TSIZE = SIZE until convergence or a singular 
matrix occurs. 

(2) TSIZE = XSIZE if TRACE * 0 after first con- 
vergence, or if a singular matrix or new components 
occur. 

(3) ITN = maximum number of iterations. 

Default: 

(1) SIZE = -In 10~ 8 = 18.420681 

(2) XSIZE = -in 10" 11 = 25.328436 

(3) ESIZE = -In 10“ 14 = 32.236191 

(4) ITN = 50 

(5) TRACE = 0 

Nondefault: 


can cause large corrections. The first situation occurs in 
the early stages of the calculation and is due to poor esti- 
mates. The second may occur at later stages of the calcu- 
lation when the iteration process sometimes attempts to 
make extremely large increases in moles of species that 
are present in small amounts. An example of this second 
situation is given in Zeleznik and Gordon (1962a). In 
both of these cases a control factor A. is used to restrict 
the size of the corrections to In n (j = 1,...,NG) and n } 
(j = NG + 1,...,NS) as well as to in n and In T obtained 
by solving the equations in tables 2. 1 and 2.2. 

The numerical value of X is determined by empiri- 
cal rules that experience has shown to be satisfactory. 
For T and n, corrections are limited to a factor of e® 4 
= 1.4918. For gas-phase species two different correction 
controls are calculated that depend on the magnitude of 
the mole fractions. The logarithm of each mole fraction 
is compared with a parameter called SIZE whose default 
value is -In 10 -8 = 18.420681. If In (rij/n) > -SIZE, 
corrections to n } are limited to a factor of e 2 = 7.3891. 
For these limitations on corrections to T, n, and rij In the 
value of a control factor X] may be calculated as 

X, = - (3.1) 

max(5|Aln T\, 5|Aln n[, |Aln n^|) 

For those gaseous species for which In (nj/n) 
< -SIZE and Ain nj > 0, a control factor X 2 is defined as 


(1) If TRACE * 0, ITN = 50 + NS/2 

(2) If TRACE < 10 -8 , XSIZE = -In TRACE and 
ESIZE = -In (TRACExlO -3 ) = XSIZE + 6.9077553 

(3) If BRATIO < 10 -5 , SIZE = In 1000/BRATIO 
and XSIZE = SIZE + In 1000 (6.9077553) 

(4) If singular matrix, XSIZE = TSIZE = 80. 


The use of ESIZE to control the size of ionized 
species permitted to be present during iteration is 
discussed in the section 3.7. 


X 2 = min 


-ln^ - 9.2103404 
n 

Ain rij - Ain n 


(3.2) 


This prevents a gaseous species with a small mole 
fraction from increasing to a mole fraction greater than 
KT 4 . The control factor X to be used in equations (3.4) 
is defined in terms of X^ and A^ as 


3.3 Convergence 


X = min(l, X x , (3.3) 


The problem of convergence is discussed in 
Zeleznik and Gordon (1962a) and Gordon and McBride 
(1976). Zeleznik and Gordon (1962a) points out that the 
iteration equations sometimes give large corrections that, 
if used directly, could lead to divergence. Two situations 


A value for X is determined for each iteration. 
Whenever current estimates of composition and/or tem- 
perature are far from their equilibrium values, X will be 
less than 1. Whenever they are close to their equilibrium 
values, A. will equal 1. New estimates for composition 
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and temperature are then obtained from the correction 
equations 


In n/ +I) = In n® + A (i) (Aln n/ n (j = 1,...,NG) 
n/ +,) = n® + A®(Ai*p® (/' = NG + 1,...,NS) 
In n (i+I) = In n (i) + A (i) (Aln n) (i) 


In T® +1) = In T® + A (0 (Aln 7) (i) 


(3.4) 


where the superscript i represents the ith estimate. 

The iteration procedure is continued until correc- 
tions to composition satisfy the following criteria: 


n.lAln n.| , 

i- <; 0.5xl0‘ 5 (j = 1,...,NG) 

NS 


J An,. | 

- — >- <l 0.5xl0~ 5 (j = NG + 1,...,NS) 

NS 

w l Aln n l <; 0.5xl0" 5 

NS 

E», (3 - 5) 


For those chemical elements for which b ^ > 
l.OxKT 6 , the convergence test for mass balance is 


NS 

*o - E a a n j 


x 1.0x10 6 

(3.6a) 


(i - !,...,«) 


where the subscript “max” refers to the chemical element 
i with the largest value of b^. When temperature is a 
variable, the convergence test for temperature is 


|Aln 71 s l.OxlO 4 ( 36b) 

For a constant-entropy problem (sp, sv, or rocket), the 
following convergence test on entropy is also required: 



s 0.5x10 4 


(3.6c) 


When TRACE * 0, an additional test is used: 


71 
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(k+l) 

^ i 


^ 0.001 


(i - l,...,t) ( 3 6d ) 


where the superscript refers to the £th iteration. The 
convergence tests in equations (3.5) and (3.6) ensure 
accuracy to five places in composition when expressed 
as mole fractions. 

As pointed out in section 3.2 the maximum number 
of iterations permitted by the CEA program is 
50 + NS/2. For most of the hundreds of different kinds 
of problems that have been solved by the program, 
convergence has been obtained in fewer them this 
number. For the first point, which starts with arbitrary 
initial estimates, a typical number of iterations is 8 to 20. 
For succeeding points, which use compositions of a 
previously calculated point for initial estimates, a typical 
number of iterations is 3 to 10. In the sample problems 
given in part II, for example, the number of iterations for 
the first point was 9 to 15. In some special cases a 
problem may be singular, and solutions to the correction 
matrix are then unobtainable. Techniques for handling 
this situation are discussed in section 3.6. 


3.4 Tests for Condensed Phases 


For the first point in a schedule of points, unless 
INSERT records are used (see part II of this report), the 
program considers only gaseous species during the itera- 
tion to convergence. For each point after the first, the 
program uses the results of a previous point for its initial 
estimate. After every convergence the program automati- 
cally checks for the inclusion or elimination of con- 
densed species. 

The test is based on the minimization of Gibbs 
energy. At equilibrium, equation (2.9) is satisfied (i.e., 
8G = 0). The requirement for a condensed species j, 
which was not previously included as a possible species, 
to now be included is that its inclusion will decrease 
Gibbs energy; that is, from equation (2.9) 


dG 

dn, 


RT 


~ E n Py < 0 (3-7) 

i-l 


where the subscript c refers to a condensed species. 
Equation (3.7) is identical to the vapor pressure test used 
in Zeleznik and Gordon (1962a) when data for the gas 
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phase as well for the condensed phase are available. The 
vapor pressure test is 


/ N 

o 

v-j 


Ji . hi 

{rt} 

c 

[RT n J 


where the subscript g refers to the gas phase of the same 
species as the solid phase referred to by the subscript c. 
Equations (3.7) and (3.8) are equivalent. The last term in 
equation (3.7) is identical for the gas and condensed 
phases of the same species. From equation (2.10) this 
term for the gas phase equals (\iJRT) g , which, using the 
definition of p in equation (2. 11), leads directly to 
equation (3.8). However, the advantage of equation (3.7) 
over (3.8) is that equation (3.8) can be used only when 
the gas phase of the species corresponding to the con- 
densed phase to be tested is present, whereas equation 
(3.7) can always be used. Using equation (3.7) eliminates 
the need for the extensive programming required in 
Zeleznik and Gordon (1962a) to accommodate 
A^C^sJ), for which gas-phase data are not available. 

At most, only one new condensed species is in- 
cluded after each convergence. In the event that several 
condensed species pass the test required by equa- 
tion (3.7), only that species giving the largest negative 
change to Gibbs energy is included as a possible species 
and convergence to a new equilibrium composition is 
obtained. This process is repeated until all condensed 
species required by equation (3.7) are included. 

If, after convergence, the concentration of a 
condensed species is negative, the species is removed 
from the list of currently considered species, and conver- 
gence to a new equilibrium composition is obtained. 

3.5 Phase Transitions and Special 
Derivatives 

The calculation method is based on the assumption 
that condensed phases are pure. Therefore, the possibility 
exists of encountering phase transition between solid and 
liquid (melting points) or between two stable solid 
phases. Such transitions constitute triple points because 
three phases of the same species coexist, one gaseous 
and two condensed. Such triple points are characterized 
by a definite vapor pressure and temperature, independ- 
ent of the relative proportions of each phase. This 
characterization is shown by the fact that the iteration 


equations of table 2.1 become singular for an assigned 
temperature and pressure and the inclusion of two 
condensed phases of the same species. At a triple point, 
for a specified system pressure, the relative amounts of 
the phases can be determined only if either the enthalpy 
or the entropy is assigned. 

The program can obtain equilibrium compositions 
containing either one or two condensed phases of a 
species with or without the corresponding gas phase. 
When temperature is assigned, no problems arise as to 
which one of the two or more condensed phases of a 
species is to be considered. However, when temperature 
is a variable, such as in combustion or rocket problems, 
several possibilities need to be considered. For example, 
if a liquid phase is being considered by the program and 
the temperature at convergence is below the melting 
point, two possibilities exist. First, the solid phase might 
be substituted for the liquid phase and a new conver- 
gence obtained, or second, both liquid and solid might 
be considered simultaneously and a new convergence 
obtained. Similar possibilites exist when convergence is 
obtained with a solid phase above the melting point. 
These possibilities and methods of treating them are 
discussed in detail in Zeleznik and Gordon (1962a); the 
discussion is not repeated herein. In brief, the following 
criteria are used by the program to determine whether to 
switch one condensed phase of a species to another or 
whether to consider both simultaneously: 

(1) Liquid present at temperature T <T m : 

(a) If T m - T > 50 K, switch solid for liquid. 

(b) If T m - T < 50 K, include solid and liquid. 

(2) Solid present at temperature T > T m : 

(a) If T - T m > 50 K, switch liquid for solid. 

(b) If T - T m < 50 K, include solid and liquid. 

Similar tests apply for two solid phases of a species. 

The unusual situation of constant temperature 
during an isentropic compression or expansion process 
can occur when two phases of the same species coexist. 
The transition temperature remains constant while one 
phase is being converted to the other. Under this 
circumstance derivatives with respect to temperature are 
not defined. Thus, several of the derivatives previously 
discussed, c p , c y9 and (3ln V/d\n 7) , cannot be obtained 
in this case. As a consequence equation (2.73) cannot be 
used to obtain y s . 

In Gordon (1970) the following expression was 
derived for the situation of constant temperature and 
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entropy: 


Ys,t 



1 

pin V \ 
\3In P) 


(3.9) 


The velocity of sound for this situation is no longer 
given by equation (2.74) but rather by 


a ST \JnRTy ST 


(3.10) 


These derivatives will be used in connection with 
discontinuities in a rocket throat (see section 6.3.4). 


3*6 Singularities 

The iteration method used in this report has 
successfully handled numerous chemical systems under 
a wide variety of thermodynamic conditions. Neverthe- 
less, special procedures are required to take care of 
several situations involving singularities that would 
otherwise cause the iteration method to fail. One such 
situation is a singularity in the coefficient matrix that is 
caused when all species that contain one of the elements 
are temporarily eliminated during the iteration process. 
Another situation resulting in a singularity occurs when 
two rows of the coefficient matrix are identical. This 
happens when the ratio of the assigned elements in these 
two rows is equal to the ratio of the stoichiometric 
coefficients of these two elements in every gaseous 
species being considered during the current iteration that 
contains both elements. That is, 

— (j = 1,...,NG) (3-11) 

% b° k 

One example for which equation (3.11) applies is 
stoichiometric hydrogen and oxygen at low temperatures 
and pressures, where the only species with w • > 10 -8 is 
H 2 0(g). Another example is stoichiometric lithium and 
fluorine at low temperatures, where the only species are 
LiF, Li 2 F 2 , and Li 3 F 3 . When either of these two 
situations occurs, the program will automatically reset all 
species currently set to zero to rtj = 10 -6 . This reset 


feature can be tried twice and will usually take care of 
these causes of singularity. 

If the restart procedure just described is not suc- 
cessful and the situation represented by equation (3.1 1) 
is the cause of the singularity, the program will attempt 
an additional procedure: (1) selecting one of the larger 
species containing both elements to be one of the new 
components; (2) reducing the number of components and 
the size of the coefficient matrix by 1 ; and (3) switching 
components for the rest of the chemical system. 

Another singularity is caused when several 
condensed species are currently being considered and 
one of these species can be formed by some linear 
combination of the others. An example of this is a 
chemical system containing iron and oxygen where the 
current iteration is simultaneously considering FeO(s), 
Fe 2 0 3 (s), and Fe 3 0 4 (s). When this occurs, the program 
automatically removes one of the species involved with 
the singularity, prints an error message with this 
information, and restarts. 


3.7 Iteration Procedure and Tests 
for Ions 

The program is capable of calculating equilibrium 
properties of plasmas (mixtures containing ionized 
species) if the plasma is considered ideal. Ideal is here 
meant to imply that no coulombic interactions are 
considered. Plasma textbooks, such as Griem (1964), 
point out that effects of coulombic interactions do need 
to be considered in plasmas (by means of the Debye- 
Huckel approximation, e.g.). However, special pro- 
gramming is needed to take care of these effects. 
Therefore, because the program does not consider these 
effects, the results of calculations when ions are 
considered will be valid only for those conditions where 
the ionic density is so small that the coulombic effects 
are unimportant. 

For ions to be considered, the charge-balance 
equation 


NG 

Ew = '° (3.12) 

is required, where a € j indicates the excess or deficiency 
of electrons in the ion relative to the neutral species. For 
example, in a mole of an ionized species, a e j = -3 for 
Ar^ and +1 for O^- To prevent difficulties in matrix 
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solutions, the program automatically removes the charge- 
balance equation when, for each ionized species being 
considered, In rij/n < -ESIZE (defined in section 3.2). 

There are situations when all the previous 
convergence tests have been passed but the ion balance 
is still incorrect. A special iteration procedure was 
developed to obtain the correct ion balance for these 
situations. It consists of obtaining a value of the 
Lagrangian multiplier for ions divided by RT , based on 
the assumption that the magnitude of the ionized species 
is small relative to the un-ionized species. The initial 
estimate for tc^ is taken to be the value in storage for the 
current point or from a previous point. The iteration 
procedure consists of the following steps: 

(1) Corrections to n e are obtained from 


(2) The test for convergence is 

|AkJ <; 0.0001 (3.14) 

(3) If this convergence test is not met, new 
estimates for the composition of ionized species are 
obtained from 

(In n } r' = (In n/ + a tj A* t ( 3 . 15 ) 

where the superscript i refers to the ith iteration. The 
previous sequence of steps is repeated until equa- 
tion (3.14) is satisfied. The CEA program allows for a 
maximum of 80 ion-balance iterations, but generally 
convergence is reached in about 2 to 10 iterations. 



Chapter 4 

Thermodynamic Data 

Thermodynamic data are included with the current 
program for reaction products and reactants. The number 
of these products and reactants changes over time as new 
ones are added to the data base. The data are selected 
from a number of sources, but the principal current 
sources are Chase et al. (1985), Cox et al. (1989), 
Gurvich et al. (1989), and Marsh et al. (1988). McBride 
et al. (1993a) documents the sources and the data for 50 
reference elements plus electron gas and deuterium, 
which are presently included in the thermodynamic data 
set. These elements are discussed in section 4.1. The data 
for the atomic gases as well as for a number of diatomic 
and polyatomic gases were calculated at NASA Lewis by 
using the PAC91 computer program described in 
McBride and Gordon (1992). 

The thermodynamic data provided with the CEA 
program are in the form of least-squares coefficients (to 
be described). These data, in formatted form, are 
processed by subroutine UTHERM and stored for further 
use in unformatted form. Subroutine UTHERM and the 
format for the coefficient data are both described in 
part II of this report. 

4.1 Assigned Enthalpies 

For each species heats of formation (and, when 
applicable, heats of transition) were combined with 
sensible heats to give assigned enthalpies H°(T). By 
definition 

H°(T) = //° (298. 15) + KH°(T) - W°(298.15)] (4.1) 

We have arbitrarily assumed H°( 298. 15) = A^//°(298. 15). 
Equation (4.1) then becomes 

H°(J) = A / «°(298.15) + [H°(T) - H° (298.15)] (4.2) 


In general, H°(T) * A f H°(T) for T * 298.15 K. 

A set of reference elements must be specified in 
order for heats of formation to be unambiguously related 
to specific reactions. Included among the species for 
which thermodynamic data are on a data file are the 
following 50 reference elements plus deuterium and 
electron gas: Ag, Al, Ar, B, Ba, Be, Bi, Br 2 , C, Ca, Cd, 
Cl 2 , Co, Cr, Cu, D 2 , e", F 2 , Fe, Ge, H 2 , He, Hg, I 2 , K, 
Kr, Li, Mg, Mn, Mo, N 2 , Na, Nb, Ne, Ni, 0 2 , P, Pb, Rb, 
S, Si, Sn, Sr, Ta, Th, Ti, U, V, W, Xe, Zn, and Zr. The 
thermodynamic data for these elements are documented 
in McBride et al. (1993a). For all reference elements 
A / //°(298.15) = //°(298.15) = 0. 

Assigned enthalpies for a number of reactants are 
in a reactant data file. For noncryogenic reactants 
assigned enthalpies (heats of formation) are given at 
298.15 K. For cryogenic liquids assigned enthalpies are 
given at their boiling points. These enthalpies are usually 
obtained by subtracting the following quantities from the 
heat of formation of the gas phase at 298.15 K: the 
sensible heat between 298.15 K and the boiling point, 
the difference in enthalpy between the ideal gas and the 
real gas at the boiling point, and the heat of vaporization 
at the boiling point. 

4.2 Least-Squares Coefficients 

For each reaction species the thermodynamic 
functions specific heat, enthalpy, and entropy as 
functions of temperature are given in the form of least- 
squares coefficients. The general form of these equations 
is as follows: 


=Y'aT q ‘ ( 4 - 3 ) 

R ^ 1 
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RT RT 


R 



(4.4) 


— *= -a, 7 2 + aj 1 In 7 + 03 + a 4 - + a s — 
RT 2 3 


(4.5) 



(4.10) 


Two sets of least-squares coefficients have been 
used in the chemical equilibrium program to represent 
thermodynamic data. The “old” form has been used for 
the last 25 years (Gordon and McBride, 1976; and 
McBride et al., 1993b). It consists of a fourth-order 
polynomial for the C°yR function plus integration 
constants for enthalpy and entropy as follows: 

= 0 , + aj * OjT 1 + a 4 T 3 + a 5 T* ( 4 - 6 ) 

R 

— =a. + —7 + — 7 2 + — 7 3 + — T 4 + — ( 4 -7) 

RT 1 2 3 4 5 7 


^-a 1 ln7- + a 1 7-*|r J ^r»*^T** 0j (48) 

The “new” form for these functions will be used to 
represent the thermodynamic data with the program 
described in part II. This form consists of seven terms 
for CyR and corresponding terms for enthalpy and 
entropy as well as the integration constants a g and a 9 as 
follows: 

^ * V 2 + <h 1 "' + a 3 + a J 

+ a s T 2 + a 6 T 3 + a, T 4 ( 4 - 9 ) 


— = -a t — -aJ 1 + a 3 In 7 + a 4 7 
R 1 2 ^ 3 4 


+ a< 


T 2 

2 




+ <I 9 


(4.11) 


For gases the temperature intervals for both the old 
and new functional forms are fixed. These intervals are 
200 to 1000 K and 1000 to 6000 K for the old form 
(i.e., the fourth-order polynomial form for C ° ). The new 
form uses these same intervals plus an additional interval 
from 6000 to 20 000 K for some gases. For the con- 
densed species each phase has its own set of coefficients. 
If possible, the old form uses the same two temperature 
intervals for condensed species as for the gases, but the 
intervals are usually limited by transition points. Further- 
more, there are two intervals only if the 1000 K common 
point is within the species temperature range. Otherwise, 
there is just one. By contrast, the new functional form 
has a flexible number of intervals in order to fit the 
selected data more accurately. 

Generally, the three functions C°JR, H°/RT, and 
S°/R are fit simultaneously. The fit is constrained to 
match the functions exactly at 7 = 298.15 K. Thus, the 
least-squares coefficients reproduce heats of formation at 
7 = 298.15 K exactly. For temperature intervals that do 
not contain 7 = 298.15 K, the fit is constrained to give 
the same functional values at the common temperature 
point for any contiguous intervals. When phase transi- 
tions occur, the fit is constrained so that the difference 
in Gibbs energy is zero between the phases at the 
transition point. 
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Chapter 5 

Thermal Transport Property Data 


The CEA program provides an option for calculat- 
ing mixture viscosities and thermal conductivities. The 
formulas for these mixture properties are given in sec- 
tion 5.2. Viscosity and thermal conductivity data for 
individual gaseous species, which are required for the 
mixture calculations, are included with the current pro- 
gram. Thermal transport properties for condensed species 
are not included inasmuch as there is no feasible method 
for calculating thermal transport properties for a multi- 
phase mixture. The thermal transport data were taken 
from Svehla (1995). 


5.2 Mixture Property Data 

Thermal conductivity and specific heat of a mixture 
each consists of two parts, the so-called “frozen” and 
“reaction” contributions. This was discussed for specific 
heat in section 2.5. Equation (2.49a) shows specific heat 
as the sum of these two parts. Analogously, thermal 
conductivity can be written as 


5.1 Data for Individual Species 

The thermal transport property data provided with 
the CEA program are in the form of least-squares coeffi- 
cients. The data for each species were fitted to the 
following form, which is also used in Gordon et al. 
(1984): 


In q 
In A. 


= A In T + — + — + D 
T t 2 


(5.1) 


where T| is viscosity and X is thermal conductivity. A 
binary interaction parameter r| - is also included for some 
pairs of species in the same form as equation (5.1). 

The coefficients in equation (5.1) were generated to 
give viscosity in units of micropoise (pP) and thermal 
conductivity in units of microwatts per centimeter-kelvin 
(pW/cm-K). The order and format of the transport data 
coefficients are given in part II of this report. 


where X^, X fr , and X re are the equilibrium, frozen, and 
reaction thermal conductivities of the mixture, 
respectively. The mixture viscosity, on the other hand, 
has only one term. 

5.2.1 Viscosity and Frozen Thermal Conductivity 


As pointed out in Gordon et al. (1984), most 
approximate mixture methods have the following form 
for mixture viscosities and frozen thermal conductivities 
(also the form used in the CEA program): 
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where NM is the number of gaseous species for thermal 
transport property calculations (NM < 50), jc ( is the mole 
fraction of species i relative to the NM gaseous species 
used for thermal transport property calculations, r) ( is the 
viscosity of species i, r| mix is the viscosity of the mix- 
ture, X' is the thermal conductivity of species i, <J>,y is the 
viscosity interaction coefficient between species i and j 
in equation (5.3) (<t>,y * 4^,), and V|/ ( y is the interaction 
coefficient between species i and j in equation (5.4) (\| 

* y,,)- 

Several forms for the interaction coefficients 4>,y and 
\\ available in the literature are compared in Gordon 
et al. (1984), and the following were used: 


where NR is the total number of chemical reactions of 
gaseous species and ^■ r H) is the heat of reaction for 
reaction i at temperature T expressed as 

NM 

(J) = E«iA° (*' = U.....NR) (5-9) 

*=1 

where are the stoichiometric coefficients written for 
the chemical reactions involving species A k as follows: 

NM 

5>^ t = 0 (i = 1,2,..., NR) (5.10) 

t=i 


4>„ 


l 




\l/4 




2M, 


M i + “j, 


1/2 


(5.5) 


and 

dr.. = <b.. 


2.41(Af. - Af.)(M. - 0.142M.) 
(M, + Mf 


(5.6) 


For some pairs of species an interaction parameter r| ■ is 
available from the literature. As discussed in section 5.1, 
these values have been least squared and are included in 
the transport property data file. When this parameter is 
available, <|>- is obtained from the following equation 
rather than from equation (5.5): 


. _ Hi 2M , 

* 1 # M i + M J 


(5.7) 


The interaction parameter r| - also appears in the expres- 
sion for the reaction contribution to thermal conductivity. 
The same values of <| )- that are obtained from equa- 
tions (5.5) or (5.7) are also used in equation (5.6). 


Equations (5.8) to (5.10) are also applicable to ionizing 
gases. In this case the heat of reaction is replaced by the 
ionization potential, and ions and electrons are species. 

Note that the stoichiometric coefficients a ijc dis- 
cussed here are defined differently from the stoichio- 
metric coefficients a- discussed in chapter 2. In the 
coefficients a ik the subscripts refer to species k in 
reaction i, and the coefficient may be positive or 
negative. In the coefficients a- the subscripts refer to 
chemical element i in species j, and the coefficient is 
always positive. 

The X ri required in equation (5.8) are found by 
solving the following set of linear equations: 

™ AH°AT) 

Y.g,K, - {i - IA...JMR) (5 11) 

y=l 


where the g tJ are given by 


NM -1 NM 


Bn = E E 

l=k*l 



(5.12) 


5.2.2 Reaction Thermal Conductivity 


and where 


The reaction contribution to thermal conductivity is 
obtained in the same manner as discussed by Butler and 
Brokaw (1957) and Brokaw (1960) and as used by 
Svehla and McBride (1973). The following equation is 
used, which is applicable when local equilibrium exists 
in a mixture of reacting or ionizing gases: 
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(5.13) 


The factor in the previous equation is a collision 
cross-section ratio that is used in Svehla and McBride 
(1973) but is not used in this report. In its place we have 
substituted the value of 1.1, which, as may be seen in 
Hirschfelder et al. (1954) (table I-N for the Lennard- 
Jones potential and table VII-E for the modified 
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Buckingham 6-exp potential), is a typical value (to two 
Figures) for this parameter. The parameter r\ kl in 
equation (5.13) was briefly discussed in section 5.2.1 and 
is available for some pairs of species with the transport 
properties of individual species. For most pairs of 
species, however, it is calculated from equation (5.7) 
where <|> w is defined, with appropriate subscripts, by 
equation (5.5). The thermal conductivity of the gas 
mixture as given by equation (5.2) is then obtained by 
adding the results of equations (5.4) and (5.8). 

5.2.3 Specific Heat for Gases Only 


An alternative method of calculating the reaction 
contribution to specific heat is given in Svehla and 
Brokaw (1966) that differs from but is equivalent to 
equation (2.49c) or the reaction part of equation (2.59). 
It is analogous to the method for obtaining the reaction 
contribution to thermal conductivity (eq. (5.8)), namely 
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(5.16) 


Some functions, such as Prandtl number, involve a 
number of other parameters (e.g., specific heat, viscosity, 
and thermal conductivity). In this event, for consistency, 
the parameters involved should be based on similar 
assumptions. Inasmuch as viscosity and thermal conduc- 
tivity are calculated for a maximum of 50 gases (with no 
condensed species permitted), specific heat, when used 
with these properties, should also be based on the same 
gases. The following equation, which is similar to 
equation (2.49a) but with different subscripts, is used to 
express specific heat for transport property calculations: 


Note that the upper case X ( in equation (5. 16) are not the 
same variables as the lower case jc*. The X i are found by 
solving the following set of linear equations: 
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where the d- are given by 

NM-l NM 

E E v«| 

f=*+l 




< X k X t j 

[ x k x t J 


(5.18) 


C pfi q C P& + 


(5.14) 


When no condensed species are present and the 
same gaseous species are included in the transport 
property calculations as in the thermodynamic property 
calculations, equations (2.49a) and (5.14) will produce 
the same numerical results. Otherwise, they may yield 
different results. 

The equation for c fr is similar to equation (2.49b) 
except for being restricted to gases and having a 
different limit for the number of gases involved; that is, 
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(5.15) 


Note that the terms in equation (5.18) are the same as 
those in equation (5.12) except that the RT/PDy group 
is not present. 

5.2.4 Prandtl Number 

Prandtl numbers have application in heat transfer 
calculations. Frozen and equilibrium Prandtl numbers 
may be calculated from previously calculated properties 
by means of 


and 
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Chapter 6 

Theoretical Rocket Performance 


Before the publication of Gordon and McBride 
(1988), the Chemical Equilibrium Calculations (CEC) 
computer program described in Gordon and McBride 
(1976) could calculate theoretical rocket performance 
only for an infinite-area combustion chamber (IAC). 
Calculation of rocket performance for a finite-area com- 
bustor (FAC), presented in Gordon and McBride (1988), 
was added as an option to the Chemical Equilibrium 
with Transport Properties (CET) program in 1988. 

Figure 6.1 presents schematic cross sections of FAC 
and IAC rocket engines. Various points at which calcula- 

Infinite-area chamber, inf 



Injector, inj Throat, t 


(a) 



Figure 6.1. — Schematics of rocket combustion chamber cross 
sections, (a) Finite-area combustion chamber, (b) infinite-area 
combustion chamber. 


tions are made in the CEA program to obtain rocket 
performance are indicated in these figures. Combustion 
and throat parameters are always calculated first auto- 
matically. For the IAC model, only one combustion 
point is calculated, namely, at infinite area (inf). How- 
ever, for the FAC model, two combustion points are 
calculated, namely, at the combustion chamber inlet (or 
equivalently at the injector face, inj) and at the combus- 
tor end, c. In addition to these two combustion points for 
the FAC, a combustion calculation for an infinite-area 
combustor, indicated in figure 6. 1 (a) by the dashed line, 
is also made. The results at this fictitious point are used 
as an aid in an iteration procedure to obtain combustor 
end conditions, as discussed in section 6.4. In addition, 
the pressure at this point is used in calculating c* (see 
section 6.2.6). Throat conditions are indicated by the 
subscript t ; other nozzle exit conditions, either subsonic 
or supersonic, are indicated by the subscript e. Nozzle 
conditions are assigned as an option and may be in the 
form of assigned area ratios, pressure ratios, or both. 

6.1 Assumptions 

The calculation of theoretical rocket performance 
involves a number of assumptions. For the same propel- 
lant and operating conditions theoretical performance can 
vary depending on which assumptions are used. For this 
report most of the assumptions are the same for both the 
IAC and FAC models. These assumptions are one- 
dimensional form of the continuity, energy, and 
momentum equations; zero velocity at the combustion 
chamber inlet; complete combustion; adiabatic combus- 
tion; isentropic expansion in the nozzle; homogeneous 
mixing; ideal-gas law; and zero temperature lags and 
zero velocity lags between condensed and gaseous 
species. The chamber in the FAC model is assumed to 
have a constant cross-sectional area. In this chamber 
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combustion is a nonisentropic, irreversible process. 
During the burning process part of the energy released is 
used to raise the entropy, and the pressure drops. 
Expansion in the nozzle is assumed to be isentropic. 

Combustion conditions are obtained with the 
assumption of chemical equilibrium of the combustion 
products. For the I AC model the CEA program provides 
the option of calculating either equilibrium or frozen 
theoretical rocket performance. Equilibrium performance 
is based on the assumption of instantaneous chemical 
equilibrium during expansion in the nozzle. Frozen per- 
formance is based on the assumption that composition 
remains frozen at the combustion composition during 
expansion. For the FAC model only equilibrium per- 
formance is permitted. 

Assuming the same velocity (either zero or other- 
wise) at the combustion chamber inlet, identical thermo- 
dynamic results are obtained for the combustion inlet 
condition for both the IAC and FAC models. 


Equation (6.2) applies only for constant-area, one- 
dimensional flow. 

6.2.2 Velocity of Flow 

The combustion chamber inlet is indicated by the 
subscripts inf for the IAC model and inj for the FAC 
model. Then using these subscripts instead of 1 and 
using e instead of 2 in equation (6.3) and assuming the 
velocity at the combustion chamber inlet to be negligible 
relative to the exit velocity result in equation (6.3) 
becoming 


u 


e 




for IAC model 

(6.5) 

for FAC model 


where h is in units of joules per kilogram and u is in 
units of meters per second. 


6.2 Parameters 


6.2.3 Force 


6.2.1 Conservation Equations 

Rocket performance, as well as other fluid dynamic 
problems in the program, is based on the following 
conservation equations, which are consistent with the 
assumptions previously discussed: 

(1) Continuity: 


From the momentum principle of fluid mechanics 
the external force on a body in a steadily flowing fluid 
is due to the change of momentum of the fluid and to 
the increase in pressure forces acting on the body. For 
rocket applications this is expressed as 

F = — - + (P e - PJA e (6.6) 

S c 


pM = Pi A i u i (61) 

(2) Momentum: 

P 2 + P2«2 = P l + Pl«l 2 (6-2) 


The conversion factor g c has been introduced to allow 
for various units. For some commonly used systems of 
units, such as the cgs system or the International System 
(Goldman and Bell, 1986), g c - 1. However, in the 
English Technical System, commonly used by engineers, 

g c = 32.1740 (lbm/lbf)(ft/s 2 ) 


(3) Energy: 



= \ 



(6.3) 


Equation (6.1) describes the condition of constant mass 
flow rate, which will be given the symbol m; that is, 


m = pAu 


(6.4) 


6.2.4 Specific Impulse 

Specific impulse is defined as force per unit mass 
flow rate. From equation (6.6) 


F = + (Pe - PM, 

m g c m 


(6.7) 
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In rocket literature specific impulse is often expressed in 
English Technical System units of pounds force per 
pound mass per second. However, for those systems of 
units previously mentioned for which g c = 1, / is both 
dimensionally and numerically equal to velocity. 

In this report when the exit pressure is equal to the 
ambient pressure, specific impulse will be given the 
symbol / sp . From equation (6.7) 


6.2.8 Coefficient of Thrust 

The coefficient of thrust is defined in terms of 
previously defined parameters 

C F = — (6.13) 

F c* 


sp 



( 6 . 8 ) 


When the ambient pressure is assumed to be zero 
(vacuum conditions), specific impulse will be given the 
symbol / vac From equations (6.7) and (6.8) 


'vac = L 


m 


(6.9) 


6.2.9 Area Ratio 

The ratio of area at any exit nozzle station to area 
at the throat is obtained from the values of equa- 
tion (6.12) at these two points 


A, _ (A/m) e 
A, (A/m) t 


(6.14) 


6.2.5 Mach Number 


Mach number is defined as the ratio of velocity of 
flow to velocity of sound: 

J(~ “ ( 6 . 10 ) 

a 

Velocity of flow is given by equation (6.5). Velocity of 
sound is given by equation (2.74) (or eq. (3.10)). 


6.2.6 Characteristic Velocity 


Characteristic velocity is given the symbol c * and 
is defined as 


PuA8c 

m 


( 6 . 11 ) 


6.2.7 Area per Unit Mass Flow Rate 

From equation (6.4) 

A _ 1 

m pu e 


( 6 . 12 ) 


The terms p and u € are obtained from equations (2.1b) 
and (6.5), respectively. 


6.3 Procedure for Obtaining 

Equilibrium Rocket Performance 
for IAC Model 

The procedure for obtaining equilibrium perform- 
ance differs somewhat for the IAC and FAC models. 
The principal difference is due to only one combustion 
point being required for the IAC model (point inf in 
fig. 6.1(b)) but two combustion points being required for 
the FAC model, namely, at the inlet and exit of the finite 
chamber (points inj and c in fig. 6.1(a)). The procedure 
for the IAC model is discussed first, inasmuch as it is 
somewhat simpler. 

For the IAC model the procedure consists of first 
determining combustion properties and then determining 
exhaust properties at the throat and at other assigned 
stations, if any, in the nozzle exit. Combustion and throat 
conditions are always obtained first automatically by the 
CEA program. To obtain other nozzle conditions (either 
subsonic or supersonic), the desired points must be 
specified as part of the input in the form of assigned area 
ratios or pressure ratios. 

For the FAC model the procedure involves first 
determining combustion properties at the combustor 
inlet. This station is also referred to as the “injector (inj) 
station.” Conditions at the end of the combustor c and at 
the throat are then both determined by means of an 
iteration loop that also includes the fictitious point inf 
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(fig. (6.1(a)), as discussed in section 6.4. The two 
combustion stations inj and c and the throat are always 
calculated first automatically by the CEA program. As is 
the case for the IAC model other nozzle exit points, if 
desired, must be specified. 


is equal to the velocity of sound. The second procedure 
is used in this report. Throat pressure is determined by 
iteration. 

The initial estimate for the pressure ratio at the 
throat is obtained from the approximate formula 


6.3.1 Combustion Conditions 

For both the FAC model at station inj and the IAC 
model at station inf, the same combustion temperature 
and equilibrium compositions are obtained by the pro- 
gram for an assigned chamber pressure and the reactant 
enthalpy (in the program HP = .TRUE.). From the 
combustion compositions and temperature the combus- 
tion entropy and other combustion properties are deter- 
mined. For the IAC model the combustion entropy s jnf 
is assumed to be constant during isentropic expansion in 
the nozzle. However, combustion in a finite-area rocket 
chamber (the FAC model) is a nonisentropic process. 
The entropy increases during the combustion process 
from s in : to s c while pressure drops from P jnj to P c . For 
the FAC model the entropy at the end of the combustion 
chamber s c is held constant during isentropic expansion 
in the nozzle. 

6.3.2 Exit Conditions 
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Equation (6.15) is found in many references on rocket 
propulsion, such as Sutton and Ross (1976), but is exact 
only when y s is constant from combustion point to 
throat. Because the value of y s is not yet known at the 
throat, the value of y s from the combustion point is used 
by the CEA program in equation (6.15). It generally 
gives an excellent initial estimate. 

Equilibrium properties for s inf and for the value of 
P f calculated from equation (6.15) are obtained as for 
any exit point. From these properties u J (using 
eq. (6.5)) and a £ (using eq. (2.74)) are calculated and 
the following test for convergence is made: 



0.4 xlO' 4 
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Exit conditions include the throat conditions and 
assigned area ratios AJA t or pressure ratios P m JP e or 
Pjnj/P^. Throat conditions are always determined 
automatically by the program. Other exit conditions, on 
the other hand, are optional and, if included, will be 
calculated after throat calculations are completed. 

For an assigned pressure ratio equilibrium compo- 
sitions and exit temperature are determined for the 
pressure P corresponding to the assigned pressure ratio 
and for the combustion entropy (j inf for IAC and s in j for 
FAC). For throat and assigned area ratios iteration 
procedures are used to determine the correct pressure 
ratios. These procedures are described in sections 6.3.3 
to 6.3.7. 

After equilibrium compositions and temperature are 
obtained for an assigned pressure ratio or area ratio, all 
the rocket parameters for that point can be determined. 

6.3.3 Throat Conditions 

Throat conditions can be determined by locating the 
pressure or pressure ratio for which the area ratio is a 
minimum or, equivalently, for which the velocity of flow 


This criterion is equivalent to ensuring that at the throat 
the Mach number is within liO^xlCT 4 . 

If the convergence test is not met, an improved 
estimate of the throat pressure ratio is obtained from the 
iteration formula 
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where the subscript k indicates the kth iteration. A 
maximum of four iterations on throat pressure is 
permitted. Usually, no more than two are required. 

6.3.4 Discontinuities at Throat 

Gordon (1970) gives a special procedure for 
obtaining throat conditions when the velocity of sound 
is discontinuous at the throat. This type of discontinuity 
may occur when a transition point, such as a melting 
point, is being calculated at the throat. The solution of 
this problem requires the following equation, which 
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permits estimating the throat pressure at the melting 
point, where the solid phase just begins to appear: 


InP, = InP + 
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where the derivative is given by equation (2.75). 

6.3.5 Assigned Subsonic or Supersonic Area Ratios 

An iteration process is used in the CEA program to 
obtain the pressure ratios corresponding to the assigned 
area ratios. In general, initial estimates for pressure ratios 
in the iteration scheme are obtained from empirical 
equations derived from plots of previously calculated 
data. However, if the current as well as the previous 
assigned supersonic area ratios are greater than 2, an 
analytic expression is used to obtain the initial estimate. 
The same analytic expression is also used in all cases to 
obtain improved estimates for the pressure ratios. 

6.3.6 Empirical Formulas for Initial Estimates of 

* v**. 

Initial estimates of pressure ratios corresponding to 
subsonic area ratios are obtained from the following 
empirical formulas: 
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When an assigned supersonic area ratio requires an 
initial estimate of pressure ratio to be obtained from an 
empirical formula, the following formulas are used: 
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In equation (6.22) the value of y s is that determined for 
throat conditions. 

6.3.7 Analytic Expression for Improved Estimates of 

The equilibrium properties obtained for the initial 
and each subsequently improved estimate of P m f/P t are 
used in equations (6.23) and (6.24) to obtain the next 
improved estimates. From table I of Gordon and 
Zeleznik (1962) the following derivative can be obtained: 
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This derivative is used in the following correction 
formula to obtain an improved estimate for P jn( /P e : 
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where the subscript k refers to the kih estimate and 
where the area ratio with no iteration subscript is the 
assigned value. The iteration procedure is continued until 


<; 0.00004 (6.25) 
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with a maximum of 10 iterations permitted. Generally, 
convergence is reached within two to four iterations. 


6.4 Procedure for Obtaining 

Equilibrium Rocket Performance 
for FAC Model 

In this report the finite-area combustor is assumed 
to have a constant cross section. For this constant-area 
combustion chamber the momentum equation (6.2) can 
be written as 


(P + P«%j - (P + P « 2 ) e (6-26) 


Substituting ml A = p« (from the continuity eq. (6.4)) into 
equation (6.26) gives 



(6.27) 
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(6.29) 




If equation (6.29) is satisfied, the calculations for 
the finite-area combustor are complete for points c and t. 
Calculations are then continued if other values of 
pressure ratio and/or area ratio have been specified in the 
input. If equation (6.29) is not satisfied, an improved 
estimate for P inf is obtained as described in section 6.4.2, 
and the procedure for the points inf, t, and c is repeated 
until equation (6.29) is satisfied. 

An iteration procedure similar to that described for 
option 1 is also used for option 2. However, the 
contraction ratio is not known for option 2. Therefore, 
the iteration procedure involves starting with an 
estimated value for AJA t as well as for P— and then 
obtaining improved estimates for both P in - and AJA V 
Not surprisingly, more iterations are required for 
option 2 than for option 1, which requires improved 
estimates for P in j only. As for option 1 iteration is 
complete when equation (6.29) is satisfied. 


6.4.1 Initial Estimates for ^inf and A J A t 


When velocity at the injector face is negligible, equa- 
tions (6.26) and (6.27) reduce to 



An iteration procedure is required to satisfy 
equation (6.28). This procedure is somewhat different for 
each of the two possible input options available for FAC. 
In option 1 the contraction ratio AJA X is assigned. In 
option 2 the mass flow rate per unit combustor area m!A c 
is assigned. The iteration procedure for option 1 is 
simpler and therefore will be described first. Four points 
shown in figure 6.1 (inj, inf, c, and t) are involved in the 
iteration procedure. Thermodynamic parameters at the 
injector face are obtained by a combustion calculation. 
This point needs to be calculated only once. Starting 
with an estimated value for P inf , calculations are then 
made for points inf, t, and c (the assigned contraction 
ratio) exactly as is done for the IAC model (for the 
infinite-area combustion point, the throat, and an 
assigned subsonic area ratio). A check is then made to 
see if equation (6.28) is satisfied to within the tolerance 


A curve given in figure 3-14 of Sutton and Ross 
(1976) relates P m JP XX y A JA t f° r m assumed value 
of y= 1.2. The following empirical equation was derived 
by fitting three selected points read from the curve: 
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Equation (6.30) is used only to obtain an initial estimate 
for ^inf 

For option 1 the assigned value of the contraction 
ratio AJA t is used in equation (6.30). For option 2 an 
initial estimate of AJA t is required. This initial estimate 
is obtained from 
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Equation (6.31) was derived by starting with the 
relationship for characteristic velocity c* = P mf A t /ih 
(eq. (6.11) with g c = 1), multiplying both sides by A c , 
and using a representative but approximate value of 
c* = 2350 m/s. When equation (6.31) is substituted into 
equation (6.30), a quadratic in P inf is obtained. The 
solution of this quadratic gives the initial estimate for 
P inf . This estimate for P mf , when used in equation (6.31), 
gives the initial estimate for AJAj. If the input value of 
m!A c is so large that equation (6.31) calculates a value 
less than 1, the program will stop the calculations and 
print out the error message INPUT VALUE OF 
MDOT/A = (value) IS TOO LARGE. GIVES CR 
ESTIMATE LESS THAN 1 . 

6.4.2 Improved Estimates for P inf and A J A t 

For option 1 an improved estimate for P- nf 
(Pjnfnew) obtained by assuming that the ratio of this 
desired P inf new value to the current value of P inf is equal 
to the ratio of the assigned value of P in j (P mj a ) to the 
current value of P in j (obtained by means of eq. (6.28)). 
This assumption leads to 

p = p ijsi£ (6.32) 

inf .new inf p 

Mnj 

Equation (6.32) often gives such an excellent improved 
estimate for P inf that it need be used only once to obtain 
convergence (eq. (6.29)). For option 2 an improved 
estimate for AJA ( is required in addition to the one for 
P inf and is obtained from 


m 



For some test cases involving hydrogen and oxygen as 
propellant, approximately four iterations involving equa- 
tion (6.33) were required for convergence (eq. (6.29)). 

6.5 Procedure for Obtaining Frozen 
Rocket Performance 

The procedure for obtaining rocket performance 
assuming that composition is frozen (infinitely slow reac- 


tion rates) during expansion is simpler than that 
assuming equilibrium composition. The reason is that 
equilibrium compositions need be determined only for 
combustion conditions. After obtaining combustion con- 
ditions in the identical way described for equilibrium 
rocket performance, the remainder of the procedure is as 
follows. 

6.5.1 Exit Conditions 

Improved estimates of the exit temperature 
corresponding to some assigned P m ffP e are obtained by 
means of the following iteration formulas: 

(>" T ,) t „ * (l» T .)„ M Aln T .) t (6 34) 

where 

(AlnT,) = 5inf -"-^ (6.35) 

w 

and where k refers to the kth estimate. The initial esti- 
mate of an exit temperature is the value of temperature 
for the preceding point. The iteration procedure is 
continued until 

|Aln T e | < 0.5x10 4 (6-36) 

The maximum number of iterations permitted by the pro- 
gram is eight, although the convergence criterion of 
equation (6.36) is generally reached in two to four 
iterations. 

Phases are also considered to be frozen. Therefore, 
the program will calculate frozen rocket performance for 
assigned schedules of pressure and/or area ratios only 
until an exit temperature is reached that is 50 K below 
the transition temperature of any condensed species 
present at the combustion point. If a calculated exit tem- 
perature is more than 50 K below the transition tempera- 
ture, this point and all subsequent points in the schedule 
are ignored by the program and do not appear in the 
output. In addition, the following message is printed: 
CALCULATIONS WERE STOPPED BECAUSE NEXT 
POINT IS MORE THAN 50 DEG BELOW TEMP 
RANGE OF A CONDENSED SPECIES. 

After an exit temperature has been determined, all 
the rocket performance parameters (section 6.2) can be 
determined. 
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6.5.2 Throat Conditions 

Calculations for frozen throat conditions are similar 
to those for equilibrium conditions. That is, equa- 
tion (6.15) is used to get initial estimates for P inf /P r 
equation (6.17) is used to get improved estimates for 
P mf /P r and equation (6.16) is used as the convergence 
criterion. With composition and phases frozen, there is 
no possibility of discontinuities at the throat, in contrast 
to equilibrium compositions. 

6.5.3 Thermodynamic Derivatives for Frozen 
Composition 

The thermodynamic derivatives discussed in pre- 
vious sections were based on the assumption that in any 
thermodynamic process, from one condition to another, 
composition reaches its equilibrium values instantan- 
eously. If, on the other hand, reaction times are assumed 
to be infinitely slow, composition remains fixed (frozen). 
In this event expressions for derivatives become simpler 
(see Zeleznik and Gordon (1968) for further discussion). 
An expression for specific heat, based on frozen compo- 


sition, has already been given in equation (2.49b). Some 
other derivatives based on frozen composition are as 
follows: From equations (2.50) and (2.51), respectively. 



(6.37) 
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From equation (2.70) 


c v = c P j ~ nR 


(6.39) 


From equations (2.73) and (6.38), it is clear that for 
frozen composition 


Y* = Y 


(6.40) 
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Chapter 7 

Incident and Reflected Shocks 


The solution to the conservation equations that 
describe conditions for incident and reflected shocks is 
most conveniently obtained in terms of assigned 
temperature before and after shock. Therefore, theoretical 
values for shock parameters, such as velocities and gas 
compositions, are generally found tabulated as functions 
of assigned temperatures (or temperature ratios). How- 
ever, in shock tube experiments shock velocities are 
generally known, rather than temperature ratios, thus 
requiring interpolation for purposes of comparison. It is 
therefore useful to have a calculating scheme that calcu- 
lates shock properties in terms of assigned velocities, and 
that is the scheme used in this report. 

No consistent nomenclature has been found in the 
literature for shock parameters. For example, in Gaydon 
and Hurle (1963) u is the velocity of the gases relative 
to the velocity of the shock front w, and v is the actual 
velocity of the gases in the shock tube. In Glass and Hall 
(1959), however, the converse is given. We will use es- 
sentially the nomenclature of Gaydon and Hurle (1963). 

Figure II. 6 in Gaydon and Hurle (1963) shows the 
usual incident and reflected conditions in both 
laboratory-fixed and shock-fixed coordinates. Velocities 
can be expressed as 


*2* = v 2 + * R (7.4) 

and 

«5 = + ** = (7.5) 

In equation (7.5) it is assumed that the shocked gases at 
the end of the shock tube are brought to rest (V 5 = 0). 
The asterisk is used with T? 2 in equation (7.4) to differ- 
entiate it from t T 2 in equation (7.3). The gas properties 
are the same for condition 2, but the relative velocities 
are different. All quantities that appear in equations (7.2) 
to (7.5) are positive, and henceforth the arrows will be 
dropped. 

7.1 Incident Shocks 

For a constant-area shock the conservation equa- 
tions describing incident conditions are those given by 
equations (6.1) to (6.3) with A j = A 2 . For iterative 
purposes it is convenient to reduce these equations to a 
form similar to that given in Zeleznik and Gordon 
(1962b) as follows: 


u=v-w = v+ w or U = w-v (7.1) 
For incident shock waves equation (7.1) gives 
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and 

£ = vv - v (7.3) 

u 2 w s y 2 

In equation (7.2) it is assumed that the test gas is at rest 
0?l = 0). For reflected waves equation (7.1) gives 
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It is also convenient to use the symbols P* and h* for 
the right sides of equations (7.6) and (7.7), respectively. 
These equations then become 
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P* = 0 


h' - hj = 0 


(7.8) 


(7.9) 
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In equation (7.9), h 2 is defined by means of 
equation (2.14). 

7.1.1 Iteration Equations 
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Applying the Newton-Raphson method to equa- 
tions (7.8) and (7.9) divided by R and using for the 
independent variables the logarithm of temperature ratio 
and pressure ratio across the shock give 
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where from equations (2.1b) and (2.3a) 


and the subscript k stands for the kth iteration. 

The partial derivatives and right sides in equa- 
tions (7.10) and (7.11) are given as follows: 
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7.1.2 Corrections and Convergence 

An arbitrary control factor X is applied to the 
corrections obtained from solution of equations (7.10) 
and (7.11) before using the corrections in equa- 
tions (7.12) to obtain improved estimates. The purpose 
of the control factor is to limit the size of the corrections 
in order to prevent divergence. Different maximum 
corrections are permitted during the course of the 
iteration process. Depending on the iteration number this 
control factor permits a maximum correction of 1.5 to 
0.05 times the previous estimates of P^P\ an d 
This is the same as permitting a maximum absolute 
correction of 0.405465 1 1 to 0.04879016 on In P*JP\ an d 
In T 2 !T v As an illustration, the control factor X, 
corresponding to the maximum correction of 0.405465 1 1 
is obtained by means of the following equations: 
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7.1.3 Initial Estimates of T^T l and PJP\ 

Formulas for temperature ratio and pressure ratio 
across the incident shock, assuming constant y, can be 
found in texts such as Gaydon and Hurle (1963). These 
formulas, slightly rearranged, are 
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A schedule of numerical values used in equations (7.20) 
and (7.21) (in addition to the value of 0.40546511) and 
the criteria for when they are used are discussed in 
part II of this report. 

Improved estimates are obtained by using equa- 
tion (7.22) with equation (7.12) to give 
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The iteration process is continued until corrections 
obtained from equations (7.10) and (7.11) meet the 
following criteria: 


If y is constant over the temperature range T ] to T 2 , 
equations (7.25) and (7.26) give exact theoretical values 
for P 2 fP\ and T 2 /T v When composition is assumed to be 
frozen across the shock, y 2 does not vary greatly from Yj, 
and equations (7.25) and (7.26) generally give excellent 
first estimates. Improved estimates are then obtained as 
described in sections 7.1.1 and 7.1.2. However, if 
equilibrium composition is assumed across the shock, 
equation (7.26) generally gives estimates of temperature 
ratio that are too high. In Gordon and McBride (1976) 
an empirical equation is used to lower the T 2 estimates 
obtained from equation (7.26). This procedure has now 
been replaced by a new estimating scheme that gives 
considerably better estimates, especially at high Mach 
numbers. The initial estimate for P 2 is still obtained from 
equation (7.25). However, the initial estimate for T 2 is 
obtained by carrying out a combustion calculation using 
the estimate for P 2 for the assigned pressure and the 
value obtained from the following equation for the 
assigned enthalpy: 

2 

K * h 2 * K + y (7-27) 

The properties for condition 1 are for the unshocked 
gases. 
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7.2 Reflected Shocks 

The reflected shock conservation equations could be 
written like equations (6.1) to (6.3), but it is more con- 
venient to use the relations in equations (7.4) and (7.5), 
which give 

Ps w * = Pz( v 2 + w *) (7.28) 

^5 + Ps W i = P 2 + P 2 ( V 2 + W R? (7.29) 

h 5 + \ w l = k 2 + \{ V 2 + w rf (7.30) 

For iterative purposes it is convenient to reduce equa- 
tions (7.28) to (7.30) to a form similar to equations (7.6) 
and (7.7) as follows: 
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For convenience, the right sides of equations (7.31) 
and (7.32) will be given the symbols F and h'. The 
reflected shock parameters may now be solved by the 
simultaneous solutions of 


P' - l* = 0 (7.33) 

P 2 

h’ - h 5 = 0 (7.34) 


The partial derivatives and right sides for equa- 
tions (7.35) and (7.36) are as follows: 




P' 


P A 


M 2 v 2 2 


Ps 

p 2 


T. 

81n — 


RT 2 /p 5 


^-1 

IP2 


(jjnTj 

^iainrj^ 


(7.40) 


7.2.1 Iteration Equations 

Applying the Newton-Raphson method to equa- 
tions (7.33) and (7.34) divided by R and using for 
independent variables the logarithm of temperature ratio 
and pressure ratio across the shock give 
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7.2.2 Corrections and Convergence 
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The same control factors, correction equations, and 
tests for convergence discussed in the previous section 
for In /y/ 5 , and In T 2 /T ] (eqs. (7.20) to (7.24)) are 
applicable for In P 5 /P 2 and In T 5 /T 2 . 
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7.2.3 Initial Estimates of T 5 fT 2 and P$!Pi 

A value of T 5 /T 2 = 2 is generally a satisfactory 
initial estimate. An estimate for P 5 /P 2 in terms of T$fT 2 
can be obtained by inverting equation (2.36) in Gaydon 
and Hurle (1963). For T 5 /T 2 = 2 this gives 
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Only one solution of equation (7.45) is positive. 



Chapter 8 

Chapman-Jouguet Detonations 


The method used for obtaining Chapman-Jouguet 
detonation parameters is described in Zeleznik and 
Gordon (1962b). There are three steps in the calculation 
procedure. First, an initial estimate of the detonation 
pressure and temperature is made. Second, the estimate 
of these parameters is improved by using a recursion 
formula. Third, the correct values are obtained by means 
of a Newton-Raphson iteration procedure. The required 
equations are derived in Zeleznik and Gordon (1962b) 
and are summarized herein for convenience in slightly 
modified form. 

The same conservation equations (6.1) to (6.3) for 
continuity, momentum, and energy that apply for shock 
also apply here with the additional constraint that 
= a 2 . For iterative purposes, as was true for the shock 
equations, it is convenient to reduce the three conserva- 
tion equations to two: 
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The pressure ratio in equation (8.1) is the reciprocal 
of the pressure ratio in equation (7.6) because, as pointed 
out in Zeleznik and Gordon (1962b), the Newton- 
Raphson iteration encounters fewer problems in the form 
of equation (8.1). 

For convenience in writing the iteration equations 
the symbols P" and h" are used to represent the right 
sides of equations (8.1) and (8.2), respectively. These 
equations then become 



h n - h 2 = 0 (8.4) 

8.1 Iteration Equations 

Applying the Newton-Raphson method to equa- 
tions (8.3) and (8.4), dividing by R, and using for 
independent variables the logarithm of temperature ratio 
and pressure ratio across the detonation give 
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where Ain PJP-, and Ain TJT, are defined by equa- 
tion (7.12). 

The partial derivatives appearing in equations (8.5) 
and (8.6) can be evaluated if y s is taken to be independ- 
ent of temperature and pressure. This assumption is 
reasonable for moderate ranges of temperatures and 
pressures. To within the accuracy of this assumption the 
partial derivatives and right sides of equations (8.5) 
and (8.6) are 


PI 


BLANK NOT HL^fTP 


PACE. 


juL 


.t^TErmCNALLY BLANK, 


39 



d( P" - — 


ain^ 


l + ,h(^y) (8.7) 

’ 2 Tw p^ainpJ rj2 


d(p” 


31n — 


P 2 (d]nV' 


Pi ^31n T)p2 


8.3 Initial Estimates of T 2 /T 1 and 

PJPl 

A good initial estimate for pressure ratio is not as 
important as a good estimate for temperature ratio. For 
a number of chemical systems that were investigated, an 
initial estimate of {P^P\\ = 15 has been found to be 
satisfactory. An initial estimate for temperature ratio is 
found by calculating the flame temperature T 2 
corresponding to the enthalpy 
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Improved initial estimates for the assumed value of 
{P^/Py ) 0 and the estimated value of (T^T ^) 0 correspond- 
ing to h 2 in equation (8.13) can be obtained by using the 
recursion formulas 
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8.2 Corrections and Convergence 

The discussion of convergence controls and correc- 
tions for shock calculations applies equally well for The quantities A/ 2 , Y 5>2 , an d c p2 in equations (8.14) 

Chapman-Jouguet detonations. Equations (7.20) to (7.22) to (8.16) are the equilibrium values for (P 2 /P {) o an( J 

give the control factor X; equation (7.23) gives the new (7' 2 /7’|) 0 . Repeating the use of equations (8.14) to (8.17) 

estimates for P^JP , and T 2 IT ] ; and equation (7.24) is the three times in the program generally provides excellent 

test for convergence. initial estimates for the Newton-Raphson iteration. 
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Chapter 9 

Input Calculations 


A number of options are provided in the program 
for specifying input, and these options are discussed in 
detail in part IT Part of the input concerns reactants. The 
total reactant may be composed of a number of reactants, 
each of which may be specified as an oxidant or a fuel. 
If the total reactant contains more than one oxidant, 
these oxidants may be combined into a total oxidant by 
specifying the relative proportion of each oxidant. 
Similarly, if the total reactant contains more than one 
fuel, the fuels may be combined into a total fuel by 
specifying the relative proportion of each fuel. The 
overall properties of the total reactant (such as elemental 
compositions b° , enthalpy h 0> internal energy Wq, 
molecular weight A/ 0 , density p 0 , positive and negative 
oxidation states V* and V~ y specific heat c 0 , and entropy 
$ 0 ) can then be calculated by specifying the relative 
amounts of total oxidant and total fuel. This method is 
particularly convenient if calculations are to be obtained 
for a number of oxidant-fuel ratios. 

To obtain assigned properties for the total reactant, 
each reactant j may be specified as either a fuel or an 
oxidant even though a reactant may ordinarily be thought 
of as inert (e.g., N 2 ). If reactants are not divided into 
fuels and oxidants, they are treated like fuels with 
o/f = 0. Letting the superscript k equal 1 for oxidant 
(provided that there is an oxidant) and 2 for fuel, the 
kilogram- atoms of the ith element per kilogram of total 
oxidant or total fuel is 

NREAC 

b, (k) = £ (i = 1 k = 1,2) (9.1) 

7=1 

where NREAC is the total number of reactants andn ; w 
is the number of kilogram-moles of the jth reactant per 
kilogram of total oxidant (k = 1) or total fuel (k = 2). If 
the amounts of oxidants and fuels are specified in terms 
of weights, rij is obtained by 
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where wj k) is the weight of the y'th reactant and Mj k> is 
the molecular weight of the jth reactant. If, however, 
amounts of oxidants and fuels are specified in terms of 

kilogram-moles N- k \ 


N m 

n, (k) = J - (j = 1,..., NREAC; k = 1,2) 

J NREAC V 

£ Nj k) Mj k) 
h i 

(9.3) 


The Mf can be calculated from the atomic weights of 
the chemical elements M i as 


K K Ea?M, 


{j = 1,..., NREAC; k = 1,2) 


(9.4) 


The bj k) can be combined by means of the 
following equation to give the kilogram-atoms per 
kilogram of total reactant: 
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where olf is the oxidant-to-fuel weight ratio. 

Formulas analogous to equations (9.1) and (9.5) can 
be used to obtain other properties of the total oxidant, 
total fuel, and total reactant (such as enthalpies, internal 
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energies, molecular weights, and densities). The energies 
given in the next four equations are in the form in which 
they appear in the matrix arrays (tables 2.1 to 2.3). 

(1) Specific enthalpy of total oxidant (k = 1) or 
total fuel ( k = 2) divided by RT, (kg-mole/kg)*^: 


(6) Molecular weight of total reactant, kg/kg-mole: 




f \ 



(9.11) 
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"RT = ^ RT ~ ^ (9 ' 6) ( If ^ = 0, M 0 = A^ 2) ; and if A/ 2) = 0, M 0 = 

(7) Density of total oxidant or fuel: 


where i s the molar enthalpy of the yth reactant, 

(J/kg-moIe)*** 

(2) Specific enthalpy of total reactant h Q divided by 
RT, kg-mole/kg: 
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(3) Specific internal energy of total oxidant 
(k = 1) or total fuel (k = 2) (u') (k) divided by RT, 
(kg-mole/kg)**^: 
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where is the molar internal energy of the yth 

reactant, (J/kg-mole)j^. 

(4) Specific internal energy of total reactant Mq 
divided by RT, kg-mole/kg: 


oo (2) 





o 

+ — 


/ 


(9.9) 


(5) Molecular weight of total oxidant (k = 1) or 
total fuel (k = 2): 
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(* = 1,2) 


(9.10) 


(8) Density of total reactant: 
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In the main program of CEA several alternative 
expressions are given for specifying the relative amounts 
of total oxidant to total fuel and for relating them to olf. 
Two expressions are given for equivalence ratio — a 
chemical equivalence ratio r and a fuel-to-air (or fuel-to- 
oxidant) equivalence ratio <|>. If the chemical equivalence 
ratio is specified, it will be necessary to make use of 

oxidation states. Let V+ and V- be positive and 
negative oxidation states of the ith element in its 
commonly occurring compounds. At least one of these 
will be zero. Thus, for example, the negative oxidation 
state for chlorine is -1 and its positive oxidation state is 
zero. The oxidation states per kilogram of total oxidant 
or total fuel are 
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The positive and negative oxidation states for the total 
reactant are then 
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The chemical equivalence ratio is now defined as 


r ■ 


V* 


V 


(9.18) 


oxidizer (e.g., H 2 0 2 ), then for off-stoichiometric 
conditions r and <(> differ for the same reaction. For 
example, for CH 3 OH + 1/2 0 2 , r = 2 and <J) = 3. 

One of the options in the SHOCK problem requires 
reactant compositions relative to the total reactant. These 
may be obtained as follows: 



for k = 1 




for Jfc = 2 


The equivalence ratio <f), commonly used in 
engineering practice, is defined as 


<t> = 



(9.19) 


U = 1,.. ..NREAC) (9 .20) 

where ntj is kilogram-moles of 7th reactant per kilogram 
of total reactant. 

Specific heat and entropy for the total reactant in 
kilogram- moles per kilogram can be calculated by using 

the Mp Cp j , and Sj of the j reactants: 


where fla is the fuel-to-air weight ratio. As Gordon 
(1982) points out, the two equivalence ratios r and <j) are 
always identical for the stoichiometric condition. In 
addition, they are always identical when all the positive 
valence atoms are in the fuel and all the negative valence 
atoms are in the oxidizer (eg., CH 4 + 0 2 ). However, if 
some negative valence atoms are in the fuel (e.g., 
CH3OH) or if some positive valence atoms are in the 
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Appendix — Symbols 


A,B,C,D 

transport property coefficients, eq. (5.1) 

A 

area, m 2 

A J A t 

ratio of nozzle exit area to throat area, 
eq. (6.14) 

A ki 

collision cross-section ratio, eq. (5.13) 

A k 

chemical formula of species k, 
dimensionless, eq. (5.10) 

a 

velocity of sound, m/s, eqs. (2.67) 
and (2.74) 

afj = 1—9) 

least-squares coefficients, eqs. (4.6) 
to (4.11) 

a ij 

stoichiometric coefficients, kilogram- 
atoms of element i per kilogram-mole of 
species j, (kg-atom)/(kg-mole) y , eq. (2.7a) 

a (t) 

y 

stoichiometric coefficients, kilogram- 
atoms of element i per kilogram-mole of 
reactant j (oxidant if k - 1, fuel if k = 2), 
(kg-atom)^V(kg-mole)j^, eq. (9.1) 

a S,T 

velocity of sound, m/s, defined by 
eq. (3.10) 


kilogram-atoms of element i per kilogram 
of mixture, (kg-atom)-/kg, eq. (2.7c) 


assigned kilogram-atoms of element i per 
kilogram of total oxidant (k = 1) or total 
fuel (k = 2), (kg-atom)fVkg (A) , eq. (9.1) 


assigned kilogram-atoms of element i per 
kilogram of total reactant, (kg-atom)/kg, 
eq. (9.5) 

c F 

coefficient of thrust, eq. (6.13) 


C°j molar heat capacity at constant pressure 

for standard state for species or reactant j, 
J/(kg-moIe) 7 (K) 

C° j molar heat capacity at constant volume 

for standard state for species j, 
J/(kg-mole) 7 (K) 

Cq specific heat at constant pressure of total 

reactant, J/kg-K, eq. (9.21) 

c p specific heat at constant pressure of 

mixture, equilibrium or frozen, J/kg-K 

c equilibrium specific heat at constant 

pressure, J/kg-K, eq. (2.49a) 

c p equilibrium specific heat at constant 

pressure for transport properties, 

J/kg-K, eq. (5.14) 

c , frozen specific heat at constant pressure, 

P ' J/kg-K, eq. (2.49b) 

c p fr frozen specific heat at constant pressure 

for transport properties, J/kg-K, eq. (5.15) 

c reaction specific heat at constant pressure, 

J/kg-K, eq. (2.49c) 

Cp re reaction specific heat at constant pressure 

for transport properties, J/kg-K, eq. (5.16) 

c v specific heat at constant volume of 

mixture, equilibrium or frozen, J/kg-K, 
eq. (2.70) or (6.39) 

c characteristic velocity, m/s, eq. (6.11) 

'j 

D kt binary diffusion coefficient, m /s, 

eq. (5.13) 

dy matrix coefficient in eq. (5.17) 
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F Helmholtz energy of mixture with 

constraints, defined by eq. (2.33), J/kg; or 
force, N, eq. (6.6) 

%F percent of total fuel in total reactant by 

weight or mass 

/ Helmholtz energy of mixture, J/kg, 

eqs. (2.30) and (2.31) 

G Gibbs energy of mixture with constraints, 

J/kg, defined by eq. (2.8) 

g Gibbs energy of mixture, J/kg, eq. (2.5) 

g c conversion factor, eqs. (6.6) to (6.8) 

HJ molar standard-state enthalpy for 

species j, J/(kg-moIe)- 

( HJ )(** molar enthalpy of reactant j (k = 1 for 

oxidant, k = 2 for fuel), J/(kg-mole)j*\ 
eq. (9.6) 

heat of formation at temperature T, 
J/kg-mole, eq. (4.2) 

A r //°(7) heat of reaction at temperature T, 
J/(kg-mole), eq. (5.9) 

h specific enthalpy of mixture, J/kg, 

eq. (2.14) 

specific enthalpy of total oxidant (k = 1) 
or total fuel ( k = 2), (J/kg) w , eq. (9.6) 

/iq specific enthalpy of total reactants, J/kg, 

eq. (9.7) 

h' term defined by right side of eq. (7.32) 

h" term defined by right side of eq. (8.2) 

h* term defined by right side of eq. (7.7) 

/ specific impulse, N/(kg/s) or m/s, 

eq. (6.7) 

/ sp specific impulse with exit and ambient 

pressures equal, N/(kg/s) or m/s, eq. (6.8) 

/ vac vacuum specific impulse, N/(kg/s) or m/s, 

eq. (6.9) 

M molecular weight of mixture, kg/kg-mole, 

eqs. (2.3) 

A/ f atomic weight of chemical element i, 

(kg/kg-atom) t , eq. (9.4) 


M, 

molecular weight of species j, eqs. (2.3b) 
and (2.4) 

M 0 

molecular weight of total reactant, 
kg/kg-mole, eq. (9.11) 

A#> 

molecular weight of total oxidant (k = 1) 
or total fuel (k = 2), (kg/kg-mole)^, eq. 
(9.10) 

A/f 

molecular weight of reactant j (oxidant if 
k = 1 ; fuel if k = 2, (kg/kg-mole)f , 
eqs. (9.4) 

MW 

molecular weight of mixture, kg/kg-mole, 
eqs. (2.4) 

J( 

Mach number, eq. (6.10) 


kilogram-moles of reactant j per kilogram 
of total reactant, (kg-mole^/kg, eq. (9.20) 

m 

mass flow rate, kg/s, eq. (6.4) 

A/f 

kilogram-moles of reactant j (oxidant if 
k = 1; fuel if k = 2), (kg-mole)f 

n 

kilogram-moles of mixture per unit mass, 
kg-mole/kg, eq. (2.2) 

n i 

kilogram-moles of species j per kilogram 
of mixture, (kg-mole)y/kg, eq. (2.2) 

nf 

kilogram-mole per kilogram of reactant j 
(oxidant if k = 1; fuel if k = 2), 
(kg-mole)f/kg rt) , eq. (9.2) or (9.3) 

o/f 

oxidant-to-fuel weight (or mass) ratio 

P 

pressure, N/m 2 

PJPe 

ratio of combustion pressure to exit 
pressure for IAC model 

P 0 

assigned or initial pressure, N/m 2 

F 

term defined by right side of eq. (7.31) 

r 

term defined by right side of eq. (8.1) 

p* 

term defined by right side of eq. (7.6) 

Pr 

Prandtl number, eqs. (5.19) 

R 

universal gas constant, 8314.51 J / 
(kg-mole)(K) 

r 

equivalence ratio based on oxidation 
states, eq. (9.18) 
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S: molar entropy of species j, 

J/(kg-mole)y(K), eq. (2.17) 

S" molar standard-state entropy for species j, 

J/(kg-mole)/K) 

s specific entropy of mixture, J/kg-K, 

eq. (2.16) 

j 0 specific entropy of total reactants, J/kg-K, 

eq. (9.22) 

T temperature, K 

UJ molar standard-state internal energy for 

species j, J/kg-mole, eq. (2.38) 

( UJ )<*) molar internal energy of reactant j 

(oxidant if k = 1; fuel if k = 2), 
(J/kg-mole)j*\ eq. (9.8) 

u velocity (in shock problems, velocity 

relative to incident or reflected shock 
front), m/s, eqs. (6.5) and (7.1) 

u' specific internal energy of mixture, J/kg, 

eq. (2.38) 

«'<*> specific internal energy of total oxidant 

( k = 1) or total fuel (k = 2), (J/kg)^, 
eq. (9.8) 

Uq specific internal energy of total reactant, 

J/kg, eq. (9.9) 

V specific volume, m /kg, eq. (2.1a) 

V* positive oxidation state of total reactant, 

eq. (9.16) 

V~ negative oxidation state of total reactant, 

eq. (9.17) 

positive oxidation state of total oxidant 
(k = 1) or total fuel ( k = 2), eq. (9.14) 

V~ ^ negative oxidation state of total oxidant 

(k = 1) or total fuel ( k = 2), eq. (9.15) 

v actual velocity of gases in shock tube, 

m/s, eq. (7.1) 

Vffi weight of jth reactant (oxidant if k - 1 ; 

fuel if k = 2), kg f\ eq. (9.2) 

w shock front velocity, m/s, eq. (7.1) 

w R reflected shock front velocity, m/s, 

eq. (7.4) 


vvy incident shock front velocity, m/s, 

eq. (7.2) 

X i unknowns in eq. (5.16) 

x- mole fraction of species i relative to all 

species (gaseous and condensed) in eq. 
(2.4a) and relative to NM gaseous species 
used for thermal transport property 
calculations in eqs. (5.3), (5.4), (5.12), 
(5.15), (5.16), and (5.18) 

a ik stoichiometric coefficients of species k in 

reaction i for transport property 
calculations, eq. (5.10) 

a k parameter defined by equation (8.16) 

y specific heat ratio, eq. (2.72) 

y v isentropic exponent, eqs. (2.71) and (2.73) 

T| i viscosity of species i, pP, eq. (5.3) 

T| •• viscosity interaction parameter, pP, 

lJ eq. (5.7) 

T) m i x viscosity of mixture, pP, eq. (5.3) 

X control factor, eqs. (3.3) and (7.22); 

thermal conductivity, eq. (5.1) 

X i Lagrangian multipliers, eqs. (2.8) and 

(2.33); thermal conductivity of species i, 
pW/cm-K, eq. (5.4) 

Xj,^ 2 , control factors, eqs. (3.1), (3.2), (7.20), 
X 3 ,A, 4 and (7.21) 

equilibrium thermal conductivity of 
mixture, pW/cm-K, eq. (5.2) 

A. fr frozen thermal conductivity of mixture, 

pW/cm-K, eq. (5.4) 

X ri variables in eqs. (5.8) and (5.11) 

X re reaction thermal conductivity of mixture, 

pW/cm-K, eq. (5.8) 

p chemical potential of species j, 

J/(kg-mole)y, eqs. (2.11) and (2.35) 

p° standard-state chemical potential of 

species j, J/(kg-mole)^, eq. (2.1 1) 

jt Lagrangian multiplier for ions divided by 

RT 
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7 i ( Lagrangian multiplier for element i 

divided by RT, eqs. (2.18), (2.19), (2.40), 
and (2.41) 

p density of mixture, kg/m 3 , eq. (2.1b) 

p 0 density of total reactant, kg/m 3 , eq. (9.13) 

p^ density of total oxidant (k = 1 ) or total 

fuel ( k = 2), (kg/m 3 ) (<:) , eq. (9.12) 

pj*^ density of reactant j (oxidant if k = 1 ; fuel 

if* = 2), (kg/m 3 )f, eq. (9.12) 

\| fjj frozen thermal conductivity interaction 

coefficient between species i and j, 
eqs. (5.4) and (5.6) 

<J> equivalence ratio based on f/a ratios, 

eq. (9.19) 

tyjj viscosity interaction coefficient between 

species i and j, eqs. (5.5) and (5.7) 

Subscripts: 

a ambient; assigned 

c end of combustion chamber; condensed 

e exit; electrons 

g gas 

inf combustor station assuming infinite 

chamber area 

inj combustor inlet or injector station for 

finite-area combustor 

k iteration k 

m melting 


P pressure 

s entropy 

stoich stoichiometric 

T temperature 

t throat 

0 total reactant; zeroth iteration 

1,2,5 stations 

Superscripts: 

° symbol for standard state; an assigned or 

initial condition 

k 1 if oxidant, 2 if fuel 

Indices: 

{ number of chemical elements (if ions are 

considered, number of chemical elements 
plus one) 

NG number of gaseous species in thermo- 

dynamic data for current chemical system 

NM number of gaseous species for transport 

property calculations 

NR number of reactions for transport property 

calculations 

NREAC number of reactants 

NS number of species, gaseous and con- 

densed, in thermodynamic data for current 
chemical system 
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TABLE 2.1— ITERATION EQUATIONS TO DETERMINE EQUILIBRIUM COMPOSITIONS FOR EITHER ASSIGNED TEMPERATURE AND PRESSURE, 
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Column not used for assigned temperature and volume. 
b Row used only for assigned enthalpy and pressure. 
c Row used only for assigned entropy and pressure. 




TABLE 2.2— ITERATION EQUATIONS TO DETERMINE EQUILIBRIUM COMPOSITIONS FOR EITHER ASSIGNED TEMPERATURE AND VOLUME, 
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‘Column not used for assigned temperature and volume. 
*Row used only for assigned internal energy and volume. 
Row used only for assigned entropy and volume. 




TABLE 2.3— EQUATIONS FOR EVALUATING DERIVATIVES WITH RESPECT TO 
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